Bundling ecosystem services at a high resolution in the UK: trade-offs and synergies in urban landscapes

Ecosystem service bundles can be defined as the spatial co-occurrence of ecosystem services in a landscape. The understanding of the delivery of multiple ecosystem services as bundles in urban areas is limited. This study modelled ecosystem services in an urban area comprising the towns of Milton Keynes, Bedford and Luton. The objectives of this study were to assess (1) how ecosystem service bundles scale at a 2 m spatial resolution and (2) identify and analyse the composition of ecosystem service bundles. Six ecosystem services were modelled with the InVEST framework at a 2 m resolution. The correlations between ecosystem services were calculated using the Spearman rank correlation coefficient method. Principal Component Analysis and K-means cluster analysis were used to analyse the distributions, spatial trade-offs and synergies of multiple ecosystem services. The results showed that regulating services had the tendency to form trade-offs and synergies. There was a significant tendency for trade-offs between supporting service Habitat quality and Pollinator abundance. Four bundle types were identified which showed specialised areas with prevalent soil erosion with high levels in water supply, areas with high values in nutrient retention, areas with high levels in carbon storage and urban areas with pollinator abundance. This study demonstrates the existence of synergies and trade-offs between ecosystem services and the formation of ecosystem service bundles in urban areas. This study provides a better understanding of the interactions between services and improve the management choices in ecosystem service provision in urban and landscape planning.


Introduction
Ecosystem services are the goods and services that ecosystems provide to society and may be classified as provisioning, supporting, regulating and cultural services (Millenium Ecosystem Assessment 2005). The ecosystem service approach involves identifying, measuring, mapping or modelling the stocks and flows of different ecosystem services and the synergies and trade-offs that may occur among them as a result of different decisions. Ecosystem service (ES) bundles have been defined as ''sets of ecosystem services that repeatedly appear together across space or time'' (Raudsepp-Hearne et al. 2010). Urbanisation, agriculture and deforestation can affect the distribution of ecosystem services and lead to spatial trade-offs and synergies between services. A spatial trade-off occurs when an ecosystem service delivery declines with respect to an increase in the delivery of another. This may be due indirectly to simultaneous response to the same driver of change on multiple ecosystem services or due to direct interactions between ecosystem services. For instance, an increase in crop production and therefore nutrient runoff could result in a decrease in water quality. Synergies occur as a result of positive interactions between ecosystem services or multiple ecosystem services can coexist. For example, restoring riparian wetlands enhances flood protection and nutrient retention which increases crop production (Bennett et al. 2009). At a landscape scale, a pattern of trade-offs between provisioning ecosystem services and both regulating and cultural ecosystem services has been found. Synergies have been found between almost all regulating services (Kong et al. 2018). It has been found that either landscape heterogeneity or deliberate management in homogeneous landscapes can enhance the delivery of ecosystem services (Crouzat et al. 2015).
Spatial trade-offs and synergies may be the result of interactions between ecosystem services due to the impact of drivers such as changing land use and climate change or the result of natural ecological processes (Fu et al. 2015). A better understanding of the relationships between ecosystem services and the affecting drivers can help in managing ecosystem services (Bennett et al. 2009). Attempts to enhance the production of one ecosystem service can result in the decline in the provision of other ecosystem services.
A number of studies have conducted assessments of ecosystem service bundles at different scales (Turner et al. 2014;Yang et al. 2015;Raudsepp-Hearne and Peterson 2016) and found that ecosystem services exhibit spatial clustering. However, ecosystem services are sensitive to scale (Grêt-Regamey et al. 2014) and to the ecosystem services selected. Studies of ecosystem service bundle analysis have been conducted at a 1, 10, and 75 km grid resolution (Turner et al. 2014;Crouzat et al. 2015;Raudsepp-Hearne and Peterson 2016). However, there are fewer studies of ecosystem service flows at a local and regional scale (de Groot et al. 2010). The aim and novelty of this study was to analyse how ecosystem service bundles are represented at a 2 m spatial resolution. Ecosystem services are produced at different scales, depend on social-ecological processes and appear to be context dependent (Andersson et al. 2015). Ecosystem services are influenced by land use and ecosystem service bundles are produced as a result of social-ecological systems (Hamann et al. 2015). Spatial bundle analysis may reveal how ecosystem services are related to human modified areas. This knowledge can assist in urban and landscape planning. Such studies address services such as climate regulation, runoff mitigation, air purification, pollination and seed dispersal, urban temperature regulation and noise reduction.
Urban green space provides multiple ecosystem services to urban areas with climate regulation, air purification, run off mitigation, cooling (Bolund and Hunhammar 1999;Gómez-Baggethun and Barton 2013) and it is deemed that the provision of ecosystem service bundles depend on the composition and configuration of urban green space (Andersson et al. 2015;Derkzen et al. 2015). The integration of ecosystem services in spatial planning shows there is a potential for green space design aimed at increasing ecosystem services (Andersson et al. 2014;Holt et al. 2015). Ecosystem services in urban areas could show different spatial patterns than in natural ecosystems as urban ecosystems are not consistent with that of natural ecosystems. Urbanisation affects the provision of ecosystem services and ecosystem functions. The influence of landscape configuration, geomorphology and hydro-ecological conditions in urban areas affects the provision of ecosystem services of nutrient cycling, soil erosion, hydrological flow causing changes in water and sediment fluxes (Alberti 2005). In this study, the Integrated Valuation of Environmental Services and Trade-offs (InVEST) modelling framework developed by the Natural Capital project was used to model six ecosystem services at a 2 m spatial resolution. The services were modelled across three urban areas: the towns of Bedford, Luton and Milton Keynes. The three towns were chosen as they exhibit a high diversity of urban forms. The combined analysis of ESs and urban structure allows the results to be more general. This enables relationships to be drawn between ESs and urban structure that could inform sustainable urban planning.
The objectives of the study were to test: (1) how ecosystem service bundles scale at a 2 m spatial resolution and (2) identify and analyse the composition of ecosystem service bundles in urban areas. This was done by modelling six ecosystem services and determining the ES gradients along which the entire bundle of the ecosystem services changed. A K-means clustering was done to identify the ES bundles. We hypothesize that by answering these objectives (1) ecosystem services will cluster to form distinct bundle types, illustrating types of spatial trade-offs and synergies between ecosystem services which when identified will improve and inform management choices and that (2) urban structure exerts an influence on the provision of ecosystem services. The results may provide insights to mitigate the spatial trade-offs and enhance synergies in the provision of ecosystem services in development planning and green space design.

Study area
The study area for this project was the combined urban areas of three large towns: Bedford, Luton and Milton Keynes, UK (Fig. 1). These towns exhibit a broad range of urban forms and histories, including historic urban centres, areas of industrial expansion and planned new town development. The focus of this study is on ecosystem services. Bedford (52°80 N, 0°270 W) originated as a medieval market town and is built on the River Great Ouse and exhibits a radial road pattern around the town centre. Its 2011 population was 106,940 and the town covers 36 km 2 , with a population density of 2971 inhabitants km -2 (Office for National Statistics 2013). Luton is a larger industrial town typified by extensive industrial parks and nineteenth century residential 'terraces' that make up much of its urban pattern (51°520 N, 0°250 W). In the 2011 census population of 258,018 and covers 58 km 2 , with a population density of 4448 inhabitants km -2 (Office for National Statistics 2013).
Milton Keynes is a planned 'new town' developed during the 1960s (52°00 N, 0°470 W), noteworthy for its unique road layout and urban form (Grafius et al. 2016). The town is structured around a grid of major roads designed for speed and ease of automotive travel, rather than the radial pattern common to many more historic English urban landscapes (Peiser and Chang 1999). The urban area possessed a population of 229,941 in 2011, covering an area of 89 km 2 with a population density of 2584 inhabitants km -2 (Office for National Statistics 2013). Milton Keynes is also characterised by a high coverage of public green space, possessing many parks and wooded foot and cycle paths (Milton Keynes Council 2015).

Land use land cover data
The fine scale (2 m) land use/land cover map used in this study was created from colour infrared aerial photography originally at 0.5 m resolution obtained from LandMap Spatial Discovery (http://landmap. mimas.ac.uk/). The imagery was taken on 2 June 2009 for Bedford, 30 June 2009 and 24 April 2010 for Luton, and 8 and 15 June 2007 and 2 June 2009 for Milton Keynes, based on cloud-free image availability. Buildings and water features were identified from UK Ordnance Survey MasterMap layers, and remaining paved surfaces were separated from vegetation through the use of a Normalised Difference Vegetation Index (NDVI) threshold. Vegetation was classified into broadleaf trees, coniferous trees, and grass/herbaceous and distinguishing vegetation types by height. Digital Terrain Models (DTMs), produced from airborne LiDAR at a 5 m resolution were used to calculate basin maps (Grafius et al. 2016). The DTMs though did not cover the study area with gaps in the data at the borders. The nutrient retention output maps therefore resulted striped at the borders. This explains the stripes at the borders on the ecosystem service bundle raster map. However, this did not affect the results and conclusions of the analysis. The gaps comprised a limited extent of the study area and a small size of the data. All raster input maps were resampled to a 2 m spatial resolution using the raster Resample tool in ArcGIS (ESRI 2017).

Ecosystem services
A total of six ecosystem service models were chosen for the towns of Bedford, Luton and Milton Keynes ( Table 1). The ecosystem services were Habitat provision, Carbon storage, Erosion control, Water provision (Quick flow and Baseflow), Nutrient retention (Nitrogen and Phosphorous retention) and Pollination. The ESs were chosen based on methods and as defined in the MA (Millenium Ecosystem Assessment 2005). The ESs are also mentioned in the main classification CICES (Haines-Young and Potschin 2018). The ESs generated by green space for the combined three urban areas allows an analysis between ecosystem service provision and the structure of the urban landscape. This approach allows the results to be more widely applicable to other urban areas across the UK. The ecosystem services were modelled with Integrated Valuation of Environmental Services and Trade-offs (InVEST) modelling suite version 3.4.4. The six models were Habitat Quality, Carbon Storage and Sequestration, Nutrient Delivery Ratio (NDR), Sediment Delivery Ratio (SDR), Sea-sonal Water Yield (SWY) and Pollination. Eight indicators were chosen as representatives of provisioning, supporting and regulating services. Carbon storage plays a role in regulating climate and has a cooling effect. In addition, carbon storage act as sinks of CO 2 . The carbon storage indicator measures the capability of carbon storage in carbon pools and the amount of carbon stored in the landscape by land cover type (A summary of the ecosystem services and indicators are shown in Table 1).

Provisioning services
The Seasonal Water Yield model was run to compute spatial indices that quantifies the relative contribution of a parcel of land to the generation of both baseflow and quick flow. Climate, soil, topography and land cover properties determine the values of baseflow and quick flow. The saturated hydraulic conductivity map  (Sharp et al. 2016). In the input table for each land use land cover (LULC) type the runoff curve number (CN) for each soil type was then assigned (USDA-NRCS 2004). The threshold flow value was set to 250,000, which coincides with a 1 km 2 contribution area (Scordo et al. 2018). The model's default parameters for alpha a parameter (0.08333), beta b parameter (1) and gamma c parameter (1) were used according to the recommendations in the InVEST user manual (Sharp et al. 2016). The baseflow maps and the quick flow maps which represent the contribution of a pixel to slow release flow and surface runoff in mm, respectively, were retained for the analysis.

Supporting services
The habitat quality model combines information on LULC suitability and threats to biodiversity to produce habitat quality maps. Habitat quality and degradation maps were derived from the habitat quality model (see Supplementary Materials).
InVEST models habitat quality as a proxy for  Baseflow is the amount of water infiltrating into the aquifer and the longer-term discharge derived from natural storages Quick flow mm Quick flow or surface runoff could give an indication of runoff mitigation by land cover type. Quick flow is the amount of surface runoff that can potentially enter the stream Notes on interpretation of ecosystem services and indicators biodiversity and for quantity and quality of available resources. A habitat to be considered of high quality must be relatively intact and function within the range of past variability. Habitat quality depends on a habitat's proximity to human land uses and the intensity of these land uses. Human land uses refer to the way in which humans employ the land and its resources. Thus, habitat quality is degraded proportionally to the rate of increase of nearby land-use. In this study for the three towns trees (broadleaf and coniferous), grass and water were considered as green space habitats and assigned a habitat score of 1 in the sensitivity table. Built-up areas were considered nonhabitat and assigned a score of 0. A habitat score is a weighting from 0 to 1 assigned to each land use/land cover type where 1 indicates the highest habitat suitability. Habitat quality scores were assigned based on biodiversity generally. The land use maps were expanded to include a buffer of 2 km using UK MasterMap data (Ordnance Survey (GB) 2017) for paved land and buildings and land cover maps (Digimap 2015). In this study, paved land and buildings were treated as 'threats'. The habitat quality model considers the effects of threats to biodiversity and habitat quality. The maximum distance over which paved land affects habitat quality was set to 2 km. The maximum distance over which buildings affect habitat quality was set to 1 km. Habitat quality was considered to be more sensitive to paved land as paved roads attract more traffic. The raster maps represent the distribution of habitat quality scores. The half-saturation constant was set to 0.5, the default for the model but can be set equal to any positive number. The half-saturation constant is a parameter and is set by the user. The InVEST model uses a half-saturation curve to convert habitat degradation scores to habitat quality scores (Sharp et al. 2016). It determines the spread and central tendency of habitat quality scores. It helps with the visual representation of heterogeneity in quality across the landscape. The habitat quality map was retained for the analysis.

Regulating services
Potential soil erosion maps were derived from the SDR model. The sediment delivery ratio model calculates erosion risks and sediment generation and flow based on topography, climate, soil and vegetation properties. The model uses a universal soil loss equation (USLE) to estimate soil erosion rates of the study area. In this study the USLE equation was adapted to UK soils. The sediment delivery model is a spatially explicit model operating at the spatial resolution of a digital elevation model (DEM) raster. For each cell, the model computes the amount of annual soil loss from that pixel, then computes the sediment delivery ratio (SDR), which is the proportion of soil loss actually reaching the stream. The model estimates soil losses and sediment delivered to the stream and retained by vegetation and topography on an annual time scale. Annual soil loss per pixel in the model is calculated according to the revised universal soil loss equation (RUSLE): where R is the rainfall erosivity (MJ mm (ha h) -1 ), K is the soil erodibility (Mg ha h (MJ ha mm) -1 ), LS is the slope length-gradient factor, C is a crop-management factor, and P is a support practise factor. The sediment delivery ratio (SDR) is calculated for each pixel i according to: where SDR i is the maximum theoretical SDR, IC is a connectivity index after Borselli et al. (2008), and IC 0 and k are calibration parameters that define the shape of the SDR function. The model generates a sediment export map which is the amount of sediment exported from each pixel that reaches the stream. Similarly, a potential soil loss map which is the amount of potential soil loss per pixel in the original land cover calculated from the USLE equation. A measure of sediment retention is produced as the difference in the amount of sediment delivered by the current land cover and a land cover where all land use types have been cleared to bare soil. For the rain erosivity (R) factor (Panagos et al. 2015), the map (REDEC ESDAC database) was resampled to a 2 by 2 m cell size. The mean values of the map were used for the study. For the soil erodobility (K) factor (Panagos et al. 2014), a measure of the susceptibility of soil particles to detachment and transport by rainfall and runoff, the map (LUCAS ESDAC database) was resampled to a 2 by 2 m grid cell size. The mean values of the map were used to run the model. For the input table the cover management factor C was parameterised for each land cover class after Morgan (2005): 0.002 for broadleaf trees; 0.004 for coniferous trees; 0.01 for grassland, 0 for buildings, water and paved surfaces. The P factor incorporates the support practices aiming at reduction of runoff from the area and refers to US farming practices. It was omitted in this study and assigned the value of 1. The threshold flow accumulation parameter was set to 250,000 (1 over the pixel area in km 2 ). The model's default parameters Borselli k parameter (2), Borselli IC 0 parameter (0.5) and maximum SDR value (0.8) were used according to the recommendations in the InVEST user manual (Sharp et al. 2016). Nitrogen (N) and phosphorous (P) retention maps were derived from the Nutrient Delivery Ratio model. The Nutrient delivery model assessed the service of nutrient retention by vegetation. The N and P retention maps were obtained using a multiplying operation between the NDR map and the nutrient loads (for surface transport) map generated by the model using the equation Nretention ¼ 1 À NDR ð ÞÂload N using the Raster Calculator in ArcGIS, respectively for nitrogen and phosphorous. The model uses a mass balance approach, describing the movement of mass of nutrient through space. The model quantifies the amounts of nitrogen and phosphorus entering watercourses with run-off and produces raster maps showing how much load from each pixel eventually reaches the stream (Sharp et al. 2016). The threshold flow accumulation was set to 250,000, which coincides with a 1 km 2 contribution area (Benez-Secanho and Dwivedi 2019). The model's default parameters for Borselli k parameter was set to 2, Subsurface critical length of Nitrogen to a global average of 200 m, Subsurface critical length of Phosphorous to a global average of 200 m, Subsurface maximum retention efficiency of Nitrogen to its default value of 0.8, the Subsurface maximum retention efficiency of Phosphorous to its default value 0.8 according to the recommendations in the InVEST user manual (Sharp et al. 2016). The nutrient retention output maps resulted striped at the borders due to missing values in the DTMs. The N and P retention maps were retained for the analysis.
Pollinator abundance was derived from the pollination model. The model was used to produce two maps, it uses estimates of the availability of nest sites and floral resources within bee flight ranges to derive an index of the abundance of bees nesting on each cell on a landscape PS (Pollinator supply) and the other showing an index of the abundance of bees visiting each cell using floral resources, and bee foraging activity and flight range information PA (Pollinator abundance). The species considered were bumble bees (Bombus spp.). For the model the input table contained information on the mean foraging distance of six species (Redhead et al. 2016). For the species B. pratorum the range foraging distance was set to 674 m (Knight et al. 2005). For each land cover a relative index of the availability of the given nesting type and floral resources was set on a scale of 0 to 1. The ground availability index for nesting suitability was set for grass to 0.5 (Grafius et al. 2018), shrub/tall grass to 1, broadleaf trees to 0.8, coniferous trees to 0.8. Paved impervious surfaces was set to 0.3 and for buildings the score was set to 0.1 for continuous urban fabric and for discontinuous urban fabric to 0.3 (Zulian et al. 2013). The pollinator abundance index maps for each of the six species were summed and divided by six using the Raster Calculator tool in ArcGIS to produce a relative mean index of pollination service provision map. The relative index mean pollinator abundance maps were retained for the analysis.
Carbon storage was calculated using the Carbon Storage and Sequestration model. The model analyses four major potential carbon sinks (1) above ground biomass (all living plant bodies); (2) root biomass (all the root systems of the above ground biomass); (3) soil organic carbon (carbon pools present in the soil); and (4) humus carbon: dead organic carbon. The model aggregates the amount of carbon stored in these pools according to the land use maps and classifications.
In the formula, C total , C above , C below , C soil and C dead represent total carbon (Mg ha -1 ), above ground biomass (Mg ha -1 ), root biomass (Mg ha -1 ), soil organic carbon (Mg ha -1 ), and humus carbon (Mg ha -1 ), respectively. A report is generated with the total amount of carbon storage in the landscape. The model generates a raster map which shows the amount of carbon currently stored in Mg in each grid cell. A study was chosen as a primary source to assign above ground biomass and the below ground biomass (Milne and Brown 1997). A study was chosen as primary source of carbon soil and humus carbon estimates of broadleaf and coniferous trees (Vanguelova et al. 2013). In the input table buildings, paved surfaces, and water were all set to zero. The total carbon storage maps for the four carbon pools were retained for the analysis.

Data analysis
JMP (SAS Institute Inc. 2018) and R (R Development Core Team 2016) software were used for statistical analysis. The ecosystem service output maps for each town were set to the extent of study area (LULC maps) and converted into tables using R. The tables were combined and the ecosystem services data were standardised by subtracting the column mean and dividing by the column standard deviation. As the ESs showed non normal distributions the correlations between ecosystem services were calculated using Spearman's rank correlation method. Subsequently a correlation test was used to assess significance levels. Negatively correlated ESs (P B 0.001) were considered trade-offs, whereas positively correlated ESs were considered synergies. Principal component analysis (PCA) was conducted in JMP to detect trade-offs and synergies among ecosystem services. The scree plot was used to determine the number of PC axes to retain by the eigenvalues that were greater than 1 (Plieninger et al. 2013) and then a factor analysis was applied to the principal components with Varimax rotation. A K-means cluster analysis was conducted in R using the 'stats' package and a maximum of 1,000 iterations was used to find the cluster solution with the lowest within-cluster sum of squares. The elbow method (Yuan and Yang 2019) was used to determine the number of clusters and the composition of ESs within the cluster were examined. The K-means function returned the mean average values and assigned cluster membership to each observation. The clusters were then converted into maps using R and used to map ecosystem service bundles for each town. The mapping of ecosystem services and bundles was done in ArcGIS Desktop 10.5.1.

Correlations between ecosystem services
Of the possible 28 pair-wise correlations, 17 pairs were significantly correlated (10 positive, 7 negative) with the strongest correlations between N retention and P retention, and N (and P) retention and carbon storage. 11 services were not correlated (r \ 0.1), 9 were weakly correlated (0.1 B r \ 0.3), 4 were moderately correlated (0.3 B r \ 0.5) and 4 were highly correlated r C 0.5). Provisioning services had primarily weak positive correlations with other services except for carbon storage (see Online Appendix, Table A1). Pollinator abundance had primarily weak and no correlations with other services and had a high negative correlation with habitat quality (r = -0.73). N retention had a high negative correlation with carbon storage (r = -0.82), in other words, places with high N retention had low values in carbon storage, and vice versa, suggesting a trade-off. P retention had a high negative correlation with carbon storage (r = -0.83). Potential soil erosion had weak positive correlations with provisioning services. Carbon storage was the ecosystem service with the highest number of negative correlations with other services and correlated negatively with six other ecosystem services.
Ecosystem service bundles PCA axes were selected according to the Kaiser-Guttman criterion, which selects the axes whose eigenvalues are greater than eigenvalue [ 1 (Legendre and Legendre, 1998). The first three PCA axes were sufficient to characterize the variation of the set of ecosystem services indicators (See Supplementary Materials for more detail, Fig. S10, S11, S12). These principal components were used to identify trade-offs and synergies among the ecosystem services. Results of the principal component analysis showed that the first three components accounted for 73.68% of the total variation (see Supplementary Materials for more detail). The first component explained 38.97% of the variance, corresponding to N retention, P retention and carbon storage. The second component explained 20.09% of the variance, corresponding to habitat quality and pollinator abundance. The third component for 14.62% of the variance, corresponding to baseflow, quick flow and potential soil erosion.
The PCA revealed trade-offs and synergies among ecosystem services as shown in the loading plot graph (see supplementary materials for more detail). Table 2 illustrates how the ecosystem services loaded among the first three PCs orthogonal axes of the component analysis factor matrix. The first axis revealed a spatial trade-off between regulating services nutrient retention and carbon storage. As values of nutrient retention increased, carbon storage decreases, indicating a trade-off. This appears to give a separation between areas with grassland with higher values in nutrient retention from areas with trees with higher values in carbon storage. The second axis showed a trade-off between habitat quality and pollinator abundance values. This appears to show a divide between suburban areas with higher values in habitat quality from urban areas with higher values in pollinator abundance, describing a gradient of stronger trade-offs from urban to suburban areas. The third factor axis represented a synergy between baseflow, quick flow and potential soil erosion as values of baseflow and quick flow increased, potential soil erosion increases. This axis is highly descriptive of baseflow and quick flow and appears to describe a gradient of stronger synergies with higher values between ESs from less built-up to more built-up areas.
The K-means cluster analysis allowed to identify four distinct bundle types for Bedford, Luton and Milton Keynes respectively (Fig. 2). The four ecosystem service bundles were clustered spatially. Bundle 1 had high above mean values in potential soil erosion and moderate values in P retention, N retention, baseflow and quick flow and comprised suburban areas where land cover is primarily grass although it also occurred in urban areas with a scattered distribution. It had above mean values for habitat quality and pollinator abundance. It had below mean values for carbon storage. It indicated areas with prevalent soil erosion ( Fig. 2 and Table 3). The key process underpinning the bundle could be water supply. Bundle 2 had high above mean values in carbon storage and moderate above mean values in habitat quality and comprised urban areas where land cover is primarily broadleaf and coniferous trees. All the other services were below the mean. It had low values for N retention and P retention, indicating areas with low values in nutrient retention (Fig. 2 and Table 3). Bundle 3 had high above mean values in pollinator abundance, nutrient retention, and comprised mainly urban areas with grassland (Fig. 3). It had above mean values for baseflow. All other services were below the mean. It had low values for habitat quality and carbon storage. Bundle 4 had the highest above mean values in nutrient retention and habitat quality with moderate above mean values for quick flow (Fig. 3). It had low values for carbon storage. It comprised suburban areas with grassland. From the cluster analysis each bundle was characterised by a higher provision of ecosystem services making it possible to name them. Bundle 1 characterised by potential soil erosion, P retention and baseflow corresponded to the (1) Potential Soil Erosion bundle; Bundle 2 characterised by carbon storage and habitat quality corresponded to the (2) Urban Trees and Woodland bundle; Bundle 3 characterised by pollinator abundance, N retention and P retention corresponded to the (3) Urban Grassland bundle; Bundle 4 characterised by N retention, P retention and habitat quality corresponded to the (4) Suburban Grassland bundle.
With regards to spatial extent the data (see Tables A2, A3 and A4) showed that for Bedford Bundle 2 had the highest spatial coverage. Bundle 3 was the most widespread with regard to areal coverage for Luton and Milton Keynes. For the combined areas Bundle 3 had the highest spatial coverage whereas Bundle 4 had the lowest spatial coverage.

Discussion
The aim of this study was to test how ecosystem service bundles scale at a 2 m resolution. It was found that ecosystem services clustered spatially to form bundles showing signs of trade-offs and synergies between ESs and distributed in specialised urban areas according to land cover type. The interactions between ESs and the synergies and trade-offs within a bundle could be explained by the linkages between ecosystem processes, functions and ecosystem services i.e.  Regulating service carbon storage had high negative correlations with nutrient retention implying a tradeoff. Such correlations most likely reflect biogeochemical processes and land cover properties. By contrast, carbon storage had a weak positive correlation only with habitat quality. Potential soil erosion had a positive correlation with baseflow and quick flow suggesting a synergy. This could be related to hydrological and climatic processes.

Bundles
The results indicate that spatial trade-offs and synergies exist among the ecosystem services indicators. The spatial trade-offs between nutrient retention on the one side and carbon storage on the other indicate a gradient of spatial trade-offs between grasslands with higher values in nutrient retention and woodland areas with higher values in carbon storage. This was the most evident separation between services revealed by the first axis of the PCA. This gradient of ecosystem service trade-offs could be explained by landscape configuration, land cover type and plant ecophysiology. Grasslands and forests have different ecological traits which affect the allocation of nutrients and carbon. Some herbaceous species are better equipped to utilize nitrate than woody species, whereas woody species are better equipped to store carbon. This is related to plant species structural and physiological characteristics and to N and C dynamics between plants, soil and climate (Groffman et al. 2004(Groffman et al. , 2009de Vries et al. 2012). In this study an evident separation between services was between suburban areas with higher values in habitat quality from urban areas with higher values in pollinator abundance revealed by the second axis in the PCA. Each bundle consists of ecosystem services that are more alike within the bundle than between bundles. The Potential Soil Erosion bundle type had the highest values for potential soil erosion and distributed in areas where land cover is primarily grass. There were signs of synergies between potential soil erosion, P retention and baseflow; a fact illustrated to some extent by the third axis in the PCA. It indicated areas with high levels of surface runoff, slow release flow, nutrient retention and potential soil erosion. This bundle was dominated by potential soil erosion. The bundle could be generated by climate and land cover properties which determine the higher levels of the ESs water supply and potential soil erosion. It had the highest number of regulating services and could indicate areas for protection (Turner et al. 2014). A similar trend was found in Chen et al. (2020) for Beijing and its surrounding areas with synergies between carbon storage and erosion control and trade-offs between erosion prevention and water supply at a coarser resolution (1 km grid cell size). A similar finding was also found by Baró et al. (2017) for the urban region of Barcelona where synergies were found between climate regulation and erosion control at 100 m resolution. Although research has found both synergic and trade-off findings for each of these ecosystem service relationships this may be related to the differences in study scales, the effects of local heterogeneity, the effects of shared ecosystem attributes (i.e. stand size, spatial configuration etc.) between ecosystem services and to the fact that urban systems are highly patchy (Wu and Li 2019).
The Urban Trees and Woodland bundle had the highest levels for carbon storage and above mean values for habitat quality. It had low values in nutrient retention and below mean values for soil erosion, quick flow, baseflow and pollinator abundance and distributed where land cover is primarily broadleaf and coniferous trees in urban areas. It indicated areas with high values in carbon storage and low values in nutrient retention. A similar trend was observed in Holt et al. (2015) for the city of Sheffield where synergies were found between carbon storage, stormwater runoff reduction and habitat provision at 500 m resolution. A comparable finding was also found in Peña et al. (2018)***** for the city of Bilbao covering an area of 413 km 2 where synergies were found between carbon storage, water flow regulation and habitat maintenance. A similar finding was found in Mexia et al. (2018) for an urban park of the city of Almada, Portugal, where synergies were found between erosion control and carbon storage at 1 m 2 resolution. A similar finding was also found in Dobbs et al. (2014) for the city of Melbourne's urban trees and woodland where synergies were found between habitat provision, climate mitigation and carbon storage at 30 m resolution. A higher supply of services was related to fragmented areas with one (or few) large patches of vegetation. These findings are consistent with previous studies that have analysed relations between ecosystem services (e.g. Derkzen et al. 2015). These studies suggest that urban green space structure plays a major role in determining the bundle of ESs.
Higher carbon storage occurs in wooded areas with trees where carbon is stored in above ground biomass (Davies et al. 2011) and soil (Edmondson et al. 2014). The interactions between ESs carbon storage and nutrient retention are most likely determined by biogeochemical processes and configuration. The low values in soil erosion and surface runoff could indicate the function of trees to anchor the soil, intercept rainfall and mitigate surface runoff and soil erosion. The Urban Grassland bundle had higher values for pollinator abundance and nutrient retention rates showing signs of spatial trade-offs with habitat quality and carbon storage and were found in urban areas where land cover is primarily grass. The low values of habitat quality are most likely due to the proximity of grassland to paved land and buildings. The Suburban Grassland bundle type had high values for nutrient retention and moderate values for habitat quality and comprised suburban areas where land cover is grassland showing signs of trade-offs with carbon storage. It indicates areas with high nutrient retention and habitat quality. The Suburban grassland bundle could be driven by land cover (grass) properties and the proximity of grassland to less built-up areas. The interactions between carbon storage and nutrient retention could be explained by plant ecological traits and functions and configuration. It appeared on the map to overlap with the Potential Soil Erosion bundle in suburban areas and had a high supply of services.
In urban ecosystems the analysis of ecosystem service bundles allows to identify the shortcomings of ecosystem service provision and has implications in land use planning (Raudsepp-Hearne et al. 2010) and allows to identify which services are vulnerable to trade-offs (Turner et al. 2014). For example, in the correlation analysis there was a weak positive correlation between habitat quality and carbon storage. However, in the Suburban grassland and Urban grassland bundle types there were signs of trade-offs between the mean values of habitat quality and carbon storage and there might be the potential to enhance carbon storage in these areas.
Previous research has found that an increase of floral resources results in an increase in pollination supply (Davis et al. 2017). In this study higher pollinator abundance was found in urban grass patches near large contiguous habitat patches of broadleaf trees and grass. This may be due to habitat suitability and proximity to likely nesting sites. Habitat connectivity enhances pollination services as pollinators benefit from high quality habitats in close proximity to nearby floral resources (Redhead et al. 2016). Allotments and residential gardens are important to pollinator species (Baldock et al. 2019). By contrast, built-up area resulted less favourable in pollinator abundance. Enhancing connectivity and floral richness is likely to support an increase in pollination services. Thus, potential for improvements in the design and implementation of environmental schemes.

Landscape structure and ecosystem service bundles
The geographical distribution and provision of ESs for the three towns appeared to be determined by urban landscape structure, land cover type and climate. The Potential soil erosion bundle occurred in suburban areas although it occurred with a scattered distribution in urban areas. The Urban trees and woodland bundle occurred in areas where land cover is trees and woodland. The Urban grassland bundle occurred in patches in urban areas where land cover is primarily grass. The Suburban grassland bundle occurred in the suburbs where land cover is primarily grass.
Moreover, the influence by land use structure and configuration on the provision of ESs was noticeable by examining the principal components maps. The PC1 map showed trade-offs between areas with high levels in nutrient retention and low levels in carbon storage (grassland) and areas with high levels in carbon storage (trees). This could be due to land use structure and configuration and plant ecological traits.
The PC2 map showed trade-offs between areas with high levels in habitat quality and low levels in pollinator abundance (suburban areas) and areas with high pollinator abundance (urban areas). For each town this appeared to separate suburban areas with higher habitat quality from urban areas with higher pollinator abundance. This appears to be related to land use with a higher variety of habitats, grassland and trees, and more connected habitats in urban areas which favours forage floral resources and nesting sites for pollinator species. For the three towns the maps showed stronger trade-offs in suburban areas than in urban areas. The increasing trade-offs are most likely due to land use and proximity to less built-up areas. The high trade-offs between habitat quality and pollinator abundance could be driven by large core areas of grassland and distance from built-up areas.
The PC3 map showed the synergies between baseflow, quick flow and soil erosion. For the three towns the synergies between ESs appeared to increase from less built-up areas to the city centre. For the town of Bedford, the map showed higher synergies between surface run-off and slow release flow in urban areas. The provision of ESs appeared to be related to urban land use structure and climate. For Luton synergies were higher in the south-east with higher values corresponding to urban built-up areas i.e. central business district. For the town of Milton Keynes, synergies between ESs were higher in urban areas and lower in the suburbs and on the east side of the town. The soil erosion and water supply synergies could be driven by climatic, topographical and hydrological processes which are likely to be influenced by landscape structure and habitat configuration. The findings may have implications in flood regulation.
The analysis of ecosystem bundles reveals differences in the amount and type of ecosystem services supplied in urban areas. These differences reflect the provision of ecosystem services in urban areas. Therefore, ES bundle analysis can be used to evaluate the supply of ESs. Landscape structure has been shown to affect the provision of multiple ecosystem services (Lamy et al. 2016). Dobbs et al. (2014) found that the provision of ESs was positively related to the amount of vegetation and negatively related to its degree of fragmentation. An urban morphology that enhances the synergies and mitigates the trade-offs of multiple ecosystem service production is important in the design of healthier and climate-resilient urban ecosystems.

Implications
ES bundle analysis can be used in urban planning to detect and evaluate the synergic and trade-off relationships between services and the variation of ES production across an urban system. It reveals how urban green space can produce bundles of ecosystem services and how land cover influences service level and is critical to green infrastructure for delivering a broad range of services for human well-being. It can provide insights to landscape planning and management, especially in the urban landscape where the interactions between patterns and processes are complex. ES bundle analysis can assist in urban planning should a decision be required to evaluate how altering urban form may influence service production. It can identify areas of high service production. For example, in this study the suburban grassland bundle had a high supply of ESs.
ES bundle analysis can be used to detect and enhance the supply of ESs. Urban green space with low levels of an ecosystem service could be enhanced. The implications for the nutrient retention service are that green space areas with high levels of trade-offs between nutrient retention and carbon storage will be represented by grasslands. The implications for the carbon storage service are that green space areas with low levels of trade-offs between nutrient retention and carbon storage will be represented by trees and woodland habitats. This gives an indication of the type of service associated with a land cover type.
The implications for the habitat quality service are that suburban green space areas with high levels of trade-offs between habitat quality and pollinator abundance will be represented mostly by grasslands and found at a distance from built-up areas. This information can give an indication of the type of service associated with a land cover type. The implications for the pollination service are that urban green space areas with low levels of habitat quality and pollinator abundance trade-offs will be found in urban built-up areas. Providing a variety of habitats in suburban areas for pollinators for nesting and forage floral resources would enhance the pollination service in those areas.
The implications for the erosion control service are that in urban green space areas with high levels of water supply and potential soil erosion synergies the erosion control service will be lower. Therefore, these areas can be closely monitored and restored. The implications for the water supply service are that green space areas with low levels in soil erosion and water supply synergies will be less at risk from surface runoff. The green space areas with high levels of soil erosion and water supply synergies suggest that these areas will have a high supply of water supply service, baseflow and quick flow.
Certain ecosystem services may be closely tied to specific land uses and associated land cover. For example, the urban trees and woodland bundle was associated with carbon storage, habitat quality and surface runoff reduction services. The landscape distribution of ecosystem services may correspond to predictable socio-ecological subsystems and associated land use types.

Limitations and uncertainties
The InVEST models were designed to be applied at a landscape scale and in non-urban environments. Therefore, there may be some bias introduced in the assessment. The sediment delivery model applies the USLE equation to calculate soil erosion which is widely used but is limited in scope, only representing rill/inter-rill erosion processes. Like the sediment delivery model, the NDR model has a small number of parameters and outputs generally show a high sensitivity to inputs. The retention efficiency values are based on empirical studies, and factors affecting these values are averaged. The model also assumes that once nutrient reaches the stream it impacts water quality at the watershed outlet, no in-stream processes are captured. The model does not represent the details of the nutrient cycle but rather represents the long-term, steady-state flow of nutrients through empirical relationships. The Seasonal Water Yield model does not include many of the complexities that occur as water moves through a landscape as reported in the manual guide. For baseflow there are no quantitative estimates only the relative contribution of pixels. Quick flow is based on curve number, which does not take topography into account. The habitat quality model assumes that all threats are additive. However, in most cases, the collective impact of multiple threats is much greater than the sum of individuals. In addition, cultural services which were not assessed could be considered in further research.

Conclusions
This study identified ecosystem service bundles and analysed the synergies and trade-offs among ecosystem services in urban areas at a 2 m spatial resolution. There was a significant tendency for trade-offs between regulating services nutrient retention and carbon storage and for regulating services to form synergies and trade-offs. There was a significant tendency for trade-offs between supporting service habitat quality and pollinator abundance. The variation in all six ESs (i.e. the entire bundle) across the landscape and the ES gradients which represent the trade-offs (or synergies) along which the ecosystem service bundle changed were revealed by the PCA. One gradient reflects areas with higher levels in nutrient retention and areas with higher levels in carbon storage. A second gradient reflects suburban areas with higher levels in habitat quality and urban areas with high levels in pollinator abundance. The six ecosystem services formed four ecosystem service bundle types showing distinct composition and geographic patterns. However, the spatial distributions of bundles are sensitive to the ecosystem services selected and care should be taken to select model data appropriate to the scale of inquiry. The bundles showed grassland areas with prevalent soil erosion and water supply, suburban grassland areas with high levels in nutrient retention and habitat quality, woodland areas with high levels in carbon storage, and urban grassland areas with high levels in pollinator abundance and nutrient retention. The results of this study show that ecosystem service bundle analysis can be used to identify the existence of synergies and trade-offs and evaluate the supply of ecosystem services. This could be used in sustainable ecosystem management and urban planning to improve the provision of multiple ecosystem services in urban areas and in green space design. Landscape composition and configuration affects the spatial distribution of ecosystem services and further studies that address the multiple relationships among ecosystem services and the influence of urban structure may provide insight and useful information to ecosystem management and to landscape planning and give a better understanding to mitigating the trade-offs and enhancing the synergies.