Introduction

Landscape ecology studies the spatial relationships between functional land units, the abiotic and biotic processes between ecosystems and the change of landscape patterns over time. Human activities are considered an integral part of ecosystems and activities such as agriculture and urban development affect the landscape structure and patterns (Mitchell et al. 2013; Qiu and Turner 2015). Landscape structure is the arrangement of land cover and land use across a landscape. Broadly, it includes landscape composition (how much of each land cover or land use that exists), configuration (the spatial pattern of these land cover or land use types), and connectivity (Mitchell et al. 2015b). Landscapes can be analysed at four different levels depending on the desired emphasis: cell, patch, class, and landscape. Landscape metrics measure landscape pattern. Patch level refers to characteristics of an individual patch, i.e. pertaining to individual patches in a categorically-classified landscape. Class level refers to a set of patches of the same type, i.e. pertaining to a single patch type (land cover type) in a categorically classified landscape. Landscape level refers to the entire patch mosaics i.e. pertaining to the full extent of the data or, as in a hierarchy, the entire patch mosaic (in a categorically classified landscape). Landscape structure or pattern is defined by its composition and spatial configuration (Leitão et al. 2006). Changes to landscape structure affect material exchange and energy flow and impede ecosystems to provide their services. Landscape structure affects the provision of multiple ESs as landscape composition and configuration affect ecological processes (Lamy et al. 2016). Urban ESs in towns and cities provide several benefits such as air purification, temperature regulation, noise reduction, urban cooling, runoff mitigation, recreational and cultural values (Bolund and Hunhammar 1999; Gómez-Baggethun and Barton 2013) which requires building a resilient supply of ESs and incorporating urban ESs in urban planning (McPhearson et al. 2015), but urban planning often entails determining trade-offs in ecosystem services and identifying where synergies occur. ES trade-offs occur when one service increases and another one decreases. Trade-offs occur when the provision of one ES is reduced as a consequence of increased use of another ES (Rodríguez et al. 2006). This may be due to simultaneous response to the same driver or due to true interactions among services. ES synergies occur when both services either increase or decrease. This may be due to simultaneous response to the same driver or due to true interactions among services (Bennett et al. 2009). In this study, the term ES interactions refer to ES trade-offs and synergies (Turkelboom et al. 2015). Modelling tools allow to map ES bundles and detect trade-offs and synergies, and within these modeling suites, Bayesian Belief Networks (BBNs) are useful tools to assess whether landscape structure affects ES trade-offs and synergies.

Significant relationships have been found between landscape structure and the provision of ecosystem services at a landscape level (Zhang and Gao 2016). Most studies have focused on landscape composition to explain the provision of ecosystem services. Lamy et al. (2016) found that landscape composition contributed more than landscape configuration in explaining the variation in the supply of ecosystem services. Incorporating landscape configuration into models though is thought to provide a better understanding of the provision of multiple ecosystem services.

Climate change and land cover patterns can affect the provision of ESs (Haase and Schwarz 2012). Urban green space and green infrastructure can enhance the resilience of cities. Urban ecosystem services include air quality regulation, climate regulation, carbon storage and sequestration, noise reduction, outdoor recreation and health benefits (Elmqvist et al. 2013; Haase et al. 2014). Dobbs et al. (2014) found synergies among regulating, provisioning and supporting services in a municipality of the city of Melbourne (in an area that covers 37.6 km2). Trade-offs were found with cultural services and regulating, provisioning and supporting services. ES provision was positively related to the amount of vegetation and negatively related to its degree of fragmentation. Peña et al. (2018) found that provisioning services had trade-offs with regulating services and cultural services in Bilbao (in an area that covers 413 km2). Synergies were found between cultural services and between regulating services. For example, synergies were found between carbon storage and water flow regulation. The trade-offs and synergies were influenced by land use type with natural land cover associated with multiple ESs. In Rotterdam, a city that covers 326 km2, more synergies between services were found than trade-offs in areas of urban green space. Synergies were found between cooling, carbon storage and air purification (Derkzen et al. 2015). The green space types (i.e. tree, woodland, shrub and herbaceous) though were not all able to provide all ESs and differed in ES provision by the amount and type of green space type in the district. ES bundles were found to depend on the composition and configuration of urban green space. Regulating services have been identified as cause for concern, and thought to underlie the sustainable production of provisioning and cultural ecosystem services and are important to the resilience of social-ecological systems (Raudsepp-Hearne et al. 2010). Verhagen et al. (2016) found that landscape configuration affects the provision of ESs at the cell and watershed scale at a 25 m resolution and that different responses of ESs to configuration suggests that accounting for configuration involves trade-offs between ESs.

Landscape structure is known to affect the provision of multiple ESs (Lamy et al. 2016) with landscape composition and configuration affecting the provision of ESs. However, there is limited knowledge on whether landscape configuration affects ES trade-offs and synergies and the mechanisms and processes that create trade-offs and synergies. Landscape configuration could drive trade-offs and synergies between ESs. The landscape configuration of urban greenspaces appears to affect their ecological and landscape functions (Leitão et al. 2006; Tian et al. 2014). Landscape configuration can affect ecological processes and the provision of ESs (Leitão et al. 2006). Landscape metrics (LM) are a useful tool to quantify landscape structure and pattern and allows to measure ecological attributes such as habitat heterogeneity, patch shape, isolation and context (Kim and Pauleit 2005), and are particularly useful in urban settings. Quantifying the spatial character of vegetation and land use can be used as a proxy for assessing the landscape’s ability to perform functions such as water and nutrient retention (Inkoom et al. 2018). Thus, the landscape structure provides an understanding of the underlying impact on ecological process.

Effective mesh size (MESH) is based on the probability that two points chosen at random will be connected. This probability is then converted into the size of a patch—the effective mesh size (unit in area, e.g. ha). The smaller the effective mesh size, the more fragmented the landscape. Effective mesh size, has been used to assess landscape fragmentation as a proxy for identifying potential habitat areas (Moser et al. 2007; Inkoom et al. 2018). Habitat fragmentation has been found to affect the provision of multiple ecosystem services (Cordingley et al. 2015) and theory predicts that fragmentation could drive trade-offs and synergies among services and may create ES bundles (Mitchell et al. 2015b). Landscape configuration in urban areas may be an important factor in influencing the provision of ESs (Holt et al. 2015). In spatial planning, LM are used in combination with other landscape pattern analytical approaches to evaluate landscape mosaics (Inkoom et al. 2018). Inkoom et al. (2018) advised that to interpret the metrics additional information about the land cover attributes and ecological functions be collected for ES assessment. Mouchet et al. (2014) provides guidelines and methods to model ES bundles and quantitative methods for analysing ES associations. Relationships between ESs can be understood by identifying which co-vary positively or negatively. Principal component analysis (PCA) can be used to identify trade-offs and synergies between services. Machine-learning algorithms, should be used to identify drivers of ES associations when the relationships among variables are complex. In this study, a BBN approach was preferred due to its ability to cope with incomplete information on the relationships between variables and uncertainty.

The aim of this study was to assess whether landscape configuration affects ecosystem service trade-offs and synergies. The objectives of this study were to test (1) a Bayesian Belief Network approach for predicting ecosystem service trade-offs and synergies in urban areas and (2) assess whether landscape configuration characteristics affect ecosystem service trade-offs and synergies. We examined the potential influence of landscape configuration on ecosystem service trade-offs and synergies by seeking answers to the following research questions: (i) does landscape configuration affect ecosystem service interactions, and (ii) can landscape configuration metrics be used to assess (and predict) ecosystem service trade-offs and synergies? We hypothesized that landscape configuration could drive ecosystem service trade-offs and synergies and that a BBN modelling approach could be used to assess the influence of landscape structure on ecosystem service trade-offs and synergies. We predict that landscape configuration affects ES trade-offs and synergies and the provision of ESs and that a BBN modelling approach can be used to test the influence of landscape configuration on ES trade-offs and synergies.

Methods

Overall approach

In this study, a BBN modelling method was used for predicting ES trade-offs and synergies and to test the relationships between landscape configuration and ES trade-offs and synergies. Patch level and class level metrics were selected to assess the influence of configuration on ES trade-offs and synergies and references based on published works suggested LM (Syrbe and Walz 2012; Haas and Ban 2018). The LM were based on a raster land use/land cover (LULC) map and were then used to assess the influence of landscape configuration on ES trade-offs and synergies by using BBN models. The sensitivity of the model outputs to the predictors was assessed in order to explore the influence of those drivers in determining predicted ES trade-offs and synergies. The ESs were modelled in the urban area comprising the towns of Milton Keynes, Bedford and Luton, UK. The input data was combined across the three towns to assess the influence of landscape configuration on ES trade-offs and synergies. This approach allows the results to be more widely applicable in other urban areas across the UK.

Study area

The study area comprised the three towns of Bedford, Luton and Milton Keynes (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. 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 km2, 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 km2, with a population density of 4448 inhabitants km−2 (Office for National Statistics 2013).

Fig. 1
figure 1

Study area showing locations and land use/land cover classification of Bedford, Luton, and Milton Keynes, UK

Full size image

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 km2 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, land structure and ecosystem services

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. Vegetated and non-vegetated surfaces were separated according to a Normalised Difference Vegetation Index (NDVI) threshold. UK Ordnance Survey MasterMap layers were used to distinguish buildings, roads and water bodies. Subsequently, airborne LiDAR was used to categorize vegetation into height classes for short grass (< 0.5 m), tall grass and shrubs (0.5–2 m), short trees (2–10 m), medium trees (10–15 m), and tall trees (> 15 m) (Grafius et al. 2016, 2019). The land cover map was resampled to a 2 m spatial resolution for all modelling and analysis. The land cover map comprised vegetation cover for broadleaf trees, coniferous trees, and grass/herbaceous and vegetation types distinguished by height (see Supplementary Materials). Note that, although the data were gathered at two different time points (2009/2010 for aerial imagery, 2012 for LIDAR), negligible change to urban land cover took place in the study area during this time.

The results of a principal component analysis (PCA) conducted in a previous study on ES datasets were used and combined with landscape metrics of mapped land cover (as described in Karimi et al. in review). The three principal components represented Nutrient retention and Carbon storage trade-offs (PC 1), Habitat quality and Pollinator abundance trade-offs (PC 2) and Potential soil erosion, Water supply synergies (PC 3). A total of six ecosystem services were modelled using the InVEST modelling framework version 3.4.4 (Sharp et al. 2016) which represented provisioning (Water supply), regulating (Carbon storage, Erosion control, Nutrient retention, Pollination) and supporting (Habitat quality) services. The PCA principal component maps (ES trade-offs and synergies) accounted for 73.68% of the total variation of the ESs. The ESs were chosen based on methods and as defined in MA (Millenium Ecosystem Assessment 2005). The ecosystem services were chosen as representatives of important provisioning, supporting and regulating services. The ecosystem services were found to cluster spatially to form four ecosystem service bundle types at a 2 m spatial resolution and exhibited distinct geographic patterns.

Landscape metrics selection and modelling method

Fragstats 4.2 (McGarigal et al. 2012) was used to calculate the following landscape metrics for the urban green spaces at the patch level; patch area (AREA), perimeter (PERIM), radius of gyration (GYRATE), perimeter-area ratio (PARA), shape index (SHAPE), fractal dimension index (FRAC), contiguity (CONTIG), core area (CORE), number of core areas in each patch (NCORE), core area index (CAI) and Euclidean Nearest-Neighbour Distance (ENN) for vegetated land cover types (see Online Appendix, Table A2). At the class level, 24 class level metrics were chosen which comprised area and edge, shape, core area, and aggregation metrics (see Table A3 for metric abbreviations and definitions). Vegetation comprised broadleaf trees, coniferous trees, and grass/herbaceous and vegetation types distinguished by height. Non-vegetated areas were treated as the background matrix and excluded from analysis. The edge depth criterion was chosen as 5 m based on research which tested core area calculations across the same study area and determined that 5 m appeared to be an effective balance across all classes (Grafius et al. 2018). For each town, the patch ID file generated by Fragstats was converted into a polygon shapefile and the principal component raster maps converted into a point shapefile. The two maps were intersected to obtain a point shapefile using the Intersect tool in ESRI ArcGIS Desktop 10.5.1. The file table was exported and joined with the calculated Fragstats patch-level metric values. The tables of each town were combined and the principal components’ values were averaged for each patch. The table was used for the analysis. The file table with the principal component map values for each town were summed and grouped by vegetation type and joined with the calculated class level metrics using JMP software (SAS Institute Inc. 2018). The tables were combined and used for the analysis.

A BBN modelling approach was chosen for its potential to use empirical data and cope effectively with incomplete information. BBNs are semi-quantitative, probabilistic models, and are useful tools for modelling ecological predictions. BBNs are suitable in an adaptive modelling framework because of their ability to update individual causal relations as new information becomes available and their explicit treatment of uncertainties in ES research. BBNs can be used to make predictions of the provision of ecosystem services and model multiple ecosystem services (Landuyt et al. 2012). An influence diagram can be used to illustrate relevance and influence between variables. It is a graph model that consists of two components: a directed acyclic graph or DAG and conditional probability tables or CPTs. The DAG consists of a set of nodes and the dependencies between the nodes are represented by directed arrows which represent cause-effect relations between dependent and independent nodes. The strengths of the causal relations between the networks variables are stored in conditional probability tables CPTs (Landuyt et al. 2013). The diagram is converted into a Bayesian Belief Network consisting of nodes and arrows showing relevance between predictor variables and response variable (Marcot et al. 2006). Furthermore, BBNs offer a range of validation techniques such as expert-based validation and sensitivity analysis (Landuyt et al. 2013). Bayesian Belief Networks have been used to model ecosystem services and in natural resource management (Cain 2001; McVittie et al. 2015).

Statistical analyses to data reduction

JMP and R software (R Development Core Team 2016) were used for statistical analysis. The reduction of landscape metrics into a limited set of metrics is necessary as redundancy between metrics can affect the outcome of the investigation. Collinearity, a high correlation between independent variables, inflates the standard errors and makes some variables less likely to be statistically significant. We used the Kolmogorov-Smirnov-Lilliefors (KSL) test and the Shapiro-Wilk test to test for data normality in JMP on the patch level and class level metrics, respectively. Further, a Spearman’s rank correlation was performed in a pair-wise correlation of 11 patch-level and 24 class-level metrics as some of our variables had a non-normal distribution using the R ‘Hmisc’ package (Harrell 2019).

Principal component factor analysis (PCA) was conducted in JMP on the patch level metrics to reduce redundancy that exists among the metrics (Hair et al. 2014). The PCA was based on the correlation matrix of the 11 metrics. The cumulative proportion of variance explained by each component was examined. An orthogonal varimax rotation was applied to redistribute the variance and facilitate the interpretation of the factor matrix. The first three principal components accounted for 76% of the total variance for all metrics (see Supplementary Materials for factor analysis results). The principal components to retain was based on the Kaiser-Guttman criterion (eigenvalue > 1). After examining the output of the correlation assessment, 7 metrics with a correlation coefficient greater or equal to \(\left| r \right|\) ≥ 0.9 were discarded by using the factor matrix derived from the PCA. The metric between the two with the lower loading was discarded. For instance, the correlation matrix revealed a high correlation between AREA and PERIM. Since AREA had a lower loading it was excluded from further analysis. Similarly, GYRATE had a high correlation with PERIM. Since GYRATE had a lower loading it was excluded. The seven patch level metrics excluded after the preliminary assessment were AREA, GYRATE, PARA, SHAPE, PERIM, NCORE and CAI. To conclude, 4 metrics were left for subsequent analysis. These included FRAC, CONTIG, CORE and ENN. To further reduce the level of redundancy between patch level metrics and identify the core patch level metrics which explained the patch variability in the dataset a second principal component factor analysis was conducted on the remaining variables and an orthogonal varimax rotation applied.

In the same way, a principal component factor analysis was conducted on the class level metrics to reduce redundancy between metrics and a varimax rotation applied. The PCA was based on the correlation matrix of the 24 metrics. The cumulative proportion of variance explained by each component was examined. The first three principal components accounted for 95% of total variance for all metrics (see Supplementary Materials). Similarly, after examining the output of the correlation assessment, 15 metrics with a correlation coefficient greater or equal to \(\left| r \right|\) ≥ 0.9 were discarded by using the factor loading matrix. The metric with the lower loading was discarded. For example, the correlation matrix revealed a high correlation between CA and PLAND. Since CA had a lower loading it was excluded. PD had a high correlation with LSI. Since PD had a lower loading it was excluded. The 15 class level metrics excluded after the preliminary assessment were CA, PD, SHAPE_AM, ED, TE, GYRATE_MN, GYRATE_AM, FRAC_AM, FRAC_MN, TCA, PLAND, CPLAND, AI, PLADJ and LPI. To conclude, 9 metrics were left for subsequent analysis. These included LSI, AREA_MN, AREA_AM, SHAPE_MN, CORE_MN, ENN_MN, ENN_AM, COHESION and MESH. To further reduce the level of redundancy among class level metrics and to identify the core landscape metrics which explained the variability in the dataset a second principal component factor analysis was conducted on the remaining variables and a varimax rotation applied (see Supplementary Materials for factor analysis results).

Model construction

BBN modelling was conducted using Netica software 6.05 (Norsys Software Corp. 2018). Separate BBN models were created with each principal component as the response variable and the landscape metrics as predictor variables. The model uses conditional probabilities to predict the response variable. Conditional probabilities define the relationships between the landscape configuration metrics and the principal component and were obtained by processing individual ‘cases’, where each case represented an ES trade-off (or synergy) value and configuration metric values found at the same location. The model then used these conditional probabilities to predict trade-offs (or synergies) at every vegetated location within the study area. All the nodes are discrete or discretised continuous variables. The case dataset with the principal component and the landscape metrics predictor variable nodes was imported and the ‘test with cases’ feature was run (Fig. 2). A sensitivity analysis was conducted on the principal component nodes. All nodes were automatically discretised with ten states each. A simplification of this with five states for the landscape metrics and three states for the response node was used for ease of visualisation of model structure and consistency when comparing conditional probabilities (Marcot et al. 2006).

Fig. 2
figure 2

Example of Bayesian Belief Network for Nutrient retention and Carbon storage trade-offs. The influence diagram illustrates the model structure and the relationship between the predictor variables, the landscape configuration a patch and b class level nodes, and the principal component node PC 1, the dependent variable. Arrows denote the direction of probabilistic influence implemented in software rather than causal relationships between the factors

Full size image

Model performance and sensitivity testing

Model performance was assessed using error rates, a goodness-of-fit measure, that expresses the frequency with which the model’s strongest prediction (most likely outcome) is incorrect against the observed data. To assess the influence of the landscape metrics in determining principal components (trade-offs or synergies between ESs) a sensitivity analysis was conducted. Sensitivity analysis determined how much the beliefs (i.e. principal component predictions) were influenced by each new finding in the predictor nodes (i.e. changes in landscape configuration characteristics). The parameter sensitivities between the predicted Nutrient retention and Carbon storage, Habitat quality and Pollinator abundance trade-offs and Potential soil erosion, Water supply synergies and the configuration metrics were assessed. Sensitivity was expressed as the expected reduction in variance of the expected real value due to a finding in a particular node (e.g. landscape metrics). The conditional probabilities for the node states were extracted and graphed as a heat map to show the predicted landscape configuration metric probability at each state level of trade-off (or synergy).

Results

Data reduction

In the Spearman correlation assessment, we found that most of the metric pairs were significant (Tables 1, 2). The first three principal components accounted for 89% of the total variance for all the patch-level metrics (Table 3). The principal components to retain was based on the percentage of variance and the Kaiser-Guttman criteria (eigenvalue > 1). Although the second and third component fell just below the threshold the first three components explained 89% of the variance and were retained as a result of percentage of variance. The scree plot shows the relationship between the increasing principal components of each metric and the cumulative proportion of variance explained (see Supplementary Materials).

Table 1 Correlation matrix showing Spearman’s rank correlation coefficients ρ (rho) illustrating the relationships between patch level landscape metrics (Significance: P < 0.001)
Full size table
Table 2 Correlation matrix showing Spearman’s rank correlation coefficients ρ (rho) illustrating the relationships between landscape class level metrics after data reduction and elimination of variables with \(\left| r \right|\) ≥ 0.9 based on factor loadings (Significance: *P < 0.05, **P < 0.01 and ***P < 0.001)
Full size table
Table 3 Result of the factor analysis for the first three factors after Varimax rotation
Full size table

As an approach to which metrics to choose, following Riitters et al. (1995) the metric with the highest loading was chosen as criterion that is representative of that factor. In summary, the metrics with the highest loadings obtained from each of our three factor loadings were: contiguity index (CONTIG, 0.89), Euclidean nearest-neighbour distance (ENN, 1), and Core area (CORE, 1), and were chosen as surrogate variables. The metrics on the first axis had the highest loadings (suggesting correlations of r > 0.7) for FRAC and CONTIG. As the metrics corresponded to shape metrics, it was termed Shape Indicator component. The metrics on the second axis had the highest loadings for ENN. As the metric corresponded to a Euclidean nearest-neighbour distance metric, it was termed Patch Distribution Indicator component. The metrics on the third axis had the highest loadings for CORE. As the metric corresponded to a core area metric, it was termed Core Area Indicator component.

Thus, among the nine remaining metrics, the first two components accounted for about 85% of the total variance of the class level metrics (Table A6). The principal components to retain was based on the Kaiser-Guttman criterion (eigenvalue > 1). As an approach to which metrics to choose, following Riitters et al. (1995) the metric with the highest loading was chosen as criterion that is representative of that factor. In summary, the metrics with the highest loadings obtained from each of our two factor loadings were: effective mesh size (MESH, 0.98) and area-weighted Euclidean nearest neighbour distance (ENN_AM, - 0.9), and were chosen as surrogate variables (Riitters et al. 1995). The scree plot shows the relationship between the increasing principal components of each metric and the cumulative proportion of variance explained. The metrics on the first axis had the highest loadings for AREA_MN, AREA_AM, CORE_MN and MESH (suggesting correlations of r > 0.7). Since the metrics corresponded to the area, core area and aggregation metrics, it was termed as Area, Core area and Aggregation Indicator component. The metrics on the second axis had the highest loadings for LSI, SHAPE_MN, ENN_MN, ENN_AM and COHESION. As the metrics corresponded to shape, shape complexity, connectivity and Euclidean nearest-neighbour distance metrics, this axis was termed Shape, Shape Complexity, Connectivity and Nearest-Neighbour Indicator component.

BBN model performance and sensitivity analysis

The LM identified in the principal component analysis were used for predicting ES trade-offs and synergies and to test the relationships between landscape configuration and ES trade-offs and synergies. At the patch level, the results of the case testing with Bayesian Belief Networks showed model error rates of 51% for PC 1, 71% for PC 2 and 67% for PC 3 (Table 4). Predicted habitat quality and pollinator abundance trade-offs and predicted potential soil erosion, water supply synergies were most sensitive to core area (CORE) and exhibited a relatively high percentage of variance reduction (respectively 20.5 and 12.8). Predicted nutrient retention and carbon storage trade-offs were most sensitive to contiguity index (CONTIG) and exhibited a high percentage in variance reduction (39.8) and had the lowest error rate (Table 4). The sensitivities were lower for ENN. This suggests that predicted nutrient retention and carbon storage trade-offs, habitat quality and pollinator abundance trade-offs and potential soil erosion, water supply synergies are least sensitive to ENN. The results indicated that CORE was the most influential predictor for two of three ES interactions at the patch level at a 2 m resolution.

Table 4 Results of case testing (error rate) and sensitivity analysis (percent in variance reduction as a metric of the relative importance of each input variable) on Bayesian Belief Network models for PC 1, PC 2 and PC 3 at a 2 m spatial resolution
Full size table

At the class level, the case testing with Bayesian Belief Networks showed an error rate of 0% for Nutrient retention and Carbon storage trade-offs, 0% for Habitat quality and Pollinator abundance trade-offs and 0% for Potential soil erosion, Water supply synergies (Table 5). The sensitivity to findings reflects the strength between predictors and principal components predictions. Predicted nutrient retention and carbon storage, habitat quality and pollinator abundance trade-offs and potential soil erosion, water supply synergies were most sensitive to MESH and exhibited a high percentage in variance reduction (respectively of 96.2, 70.7 and 97.5). This indicated that MESH was the most influential predictor in determining ecosystem service trade-offs and synergies in urban areas at the class level at a 2 m spatial resolution.

Table 5 Results of case testing (error rate) and sensitivity analysis (percent in variance reduction as a metric of the relative importance of each input variable) on Bayesian Belief Network models for PC 1, PC 2 and PC 3 at a 2 m spatial resolution
Full size table

Probabilistic associations between landscape configuration and ES trade-offs and synergies

Heat maps show the nature of probabilistic associations between the landscape configuration metrics and predicted trade-offs or synergies levels using simplified node levels (Fig. 3, Table 6). Here, high conditional probabilities reflect the likelihood of an outcome given a set of parent node states; e.g. low and moderate nutrient retention and carbon storage trade-offs are expected in areas with low core area, whereas high nutrient retention and carbon storage trade-offs are expected in patches with a large core area, a low Euclidean nearest neighbour distance and a high contiguity index. Core area (CORE) appeared to be a strong predictor for trade-offs of habitat quality and pollinator abundance, with the highest conditional probabilities associated with low core area at low and moderate levels of predicted habitat quality and pollinator abundance trade-offs. High levels of predicted potential soil erosion, water supply synergies were associated with low core area, low Euclidean nearest neighbour distance and a high contiguity index.

Fig. 3
figure 3figure 3

Heat maps that visually depict the conditional probabilities driving each model of the a patch and b class level metrics. These represent the strength of the relationships between the landscape configuration metrics and the predicted trade-offs (or synergies) of each ecosystem service trade-offs or synergies type (Bin range values are shown in Table 6). Darker cells denote higher conditional probabilities, i.e. a higher likelihood of an outcome given that set of conditions, or a stronger relationship between the combination of input parameter values and predicted trade-offs (or synergies) represented by that cell in the heat map

Full size image
Table 6 Bin ranges for input parameter values in Fig. 3 heat maps of Bayesian model conditional probabilities
Full size table

Mesh size (MESH) appeared to be a strong predictor for trade-offs of nutrient retention and carbon storage, with the highest conditional probabilities associated with low mesh size at low to moderate levels of predicted trade-offs. High levels of predicted nutrient retention and carbon storage trade-offs appeared to be associated with high mesh size. Mesh size appeared to be a strong predictor for Habitat quality and Pollinator abundance trade-offs, with the highest conditional probabilities associated with low mesh size at low to moderate levels of predicted trade-offs. High levels of predicted potential soil erosion, water supply synergies appeared to be associated with low mesh size. Low levels of predicted potential soil erosion, water supply synergies appeared to be associated with high mesh size.

Discussion

Using a BBN modelling approach this study aimed to test whether landscape configuration affects ES trade-offs and synergies. The levels in the sensitivity to findings of the dependent variables (ES interaction predictions) to the landscape metrics suggest that landscape configuration metrics core area and mesh size are the most influential determinants of ecosystem service trade-offs and synergies.

Selected landscape metrics

The first three principal components accounted for about 89% of the total amount of variance observed in the data of the patch level metrics. Out of the four metrics assessed with principal component factor analysis, factor loadings on the first axis suggested strong contributions from FRAC and CONTIG. Factor loadings on the second and third factor axes had strong contributions from ENN and CORE. AREA, GYRATE, PERIM, PARA, SHAPE, NCORE and CAI were excluded due to high correlation in the Spearman correlation assessment. Our results suggest that FRAC and CONTIG could be used to assess ecosystem services and in spatial planning. Likewise, the high loadings on the second and third factor axes suggested that ENN and CORE could be used to assess ecosystem services. However, the influence of landscape structure on ecosystem service trade-offs and synergies in urban areas can be better studied and captured using our metrics CONTIG, ENN and CORE.

FRAC and CONTIG are shape metrics. Shape metrics can be used to gain insight into the processes forming the mosaic, for example shape and shape complexity can be used for habitat services and scenic quality (Syrbe and Walz 2012). CORE has been used to assess landscape accessibility and habitat quality (Inkoom et al. 2018). ENN can be used to measure habitat accessibility. The findings suggest the possibility of using metrics CONTIG and FRAC to assess and investigate ES trade-offs and synergies in urban areas. Likewise, ENN and CORE could be used in ES assessment.

The first two components accounted for about 85% of the total amount of variance observed in the data of the class level metrics. Out of the nine metrics assessed with principal component factor analysis, factor loadings on the first axis suggested strong contributions from AREA_MN, AREA_AM, CORE_MN and MESH. CA, PD, SHAPE_AM, ED, TE, GYRATE_MN, GYRATE_AM, FRAC_AM, FRAC_MN, TCA, PLAND, CPLAND, AI, PLADJ and LPI were excluded due to high correlation in the Spearman correlation assessment. Our results suggest that AREA_MN, AREA_AM, CORE_MN and MESH could be used to assess the influence of landscape structure on ecosystem service trade-offs and synergies and in spatial planning. The metrics which loaded on the second axis LSI, SHAPE_MN, ENN_MN, ENN_AM and COHESION could be used to investigate ecosystem service trade-offs and synergies in ES assessment and urban planning. However, the influence of landscape configuration on ecosystem service trade-offs and synergies in urban areas can be better studied and captured by the use of our metrics MESH and ENN_AM.

The metrics on the second axis LSI, SHAPE_MN, ENN_MN, ENN_AM and COHESION represented shape, shape complexity, connectivity and Euclidean nearest-neighbour distance metrics. The metrics have the potential for ES assessment and land use planning. COHESION has been used to measure connectivity at the class level for pollination, pest regulation, seed dispersal, and habitat services (Haas and Ban 2018).

Predicted ES trade-offs and synergies and landscape configuration influence on ES trade-offs and synergies

Broad generalities can be established from the conditional probability heat maps about the relationships between ES trade-offs and synergies and landscape configuration (Fig. 3). At the patch level, CORE emerged as an important factor for trade-offs (nutrient retention and carbon storage, habitat quality and pollinator abundance), with low core area and low Euclidean nearest-neighbour distance associated with low and moderate levels of predicted trade-offs. The conditional probabilities show that low nutrient retention and carbon storage trade-offs predictions are expected in areas with low core area and low Euclidean distance. Low core area and low Euclidean distance could be an indication of greenspace patches characterised by loosely-connected core areas of trees, as individuals or in small stands, and vegetated road verges, wooded corridors and buffered streams. Urban woodland was associated with low nutrient retention and carbon storage trade-offs. Therefore, low core area associated with low nutrient retention and carbon storage trade-offs prediction could give an indication of fragmented woodland. The spatial distribution of the patches and relative location may matter as well (Leitão et al. 2006; Verhagen et al. 2016).

The conditional probabilities show that high habitat quality and pollinator abundance trade-offs predictions are expected in areas with a large core area. Large patches of contiguous grassland in the suburbs were associated with high levels of habitat quality and pollinator abundance trade-offs. 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 low levels of pollinator abundance could be due to low proximity to likely nesting sites for pollinator species and lower floral resources (Gunnarsson and Federsel 2014; Baldock et al. 2015).

High potential soil erosion, water supply synergies predictions are expected in areas with a low core area. The increase in the synergies could be due to a smaller patch size to mitigate surface runoff and soil erosion. Large patches of vegetation cover protect aquifers and soil resources (Kim and Pauleit 2005). This could be an effect of fragmentation and shrinkage on service provision at a landscape level. The predicted synergies could be driven by core area and by ecological processes, i.e. water flow.

At the class level MESH emerged as an important factor for nutrient retention and carbon storage trade-offs, with low MESH associated with low and moderate levels of predicted trade-offs. Mesh size can be used to measure erosion or flood prevention (Syrbe and Walz 2012; Inkoom et al. 2018). Larger values indicate a higher capacity for vegetation land cover to mitigate soil erosion and surface runoff. High surface runoff and flood risks appear to be associated with green space with a low mesh size.

In this study the urban ecosystem services were influenced by landscape structure. The findings show that urban habitat configuration exerts an influence on ES trade-offs and synergies. The surrounding matrix appeared to influence habitat configuration and ES provision. Suburban green areas were mainly associated with large patches of grassland with high nutrient retention and carbon storage trade-offs whereas urban green areas were mainly associated with smaller patches of grassland and isolated trees. The differences in trade-offs between suburban areas and urban areas could partly be driven by a higher degree of fragmentation in urban areas. The high levels of nutrient retention and carbon storage, and habitat quality and pollinator abundance predicted trade-offs could be due to land use configuration (i.e. fragmentation), land-use changes and ecological processes.

Climate change affects the distribution of ecosystems, species and processes and one of the best-documented anthropogenic climate modifications is the urban heat island (UHI) effect (Grimm et al. 2008). Climate change and land cover patterns can affect the provision of ESs (Haase and Schwarz 2012). Green infrastructure and nature-based solutions can respond to environmental change and enhance the provision of ESs and the benefits they provide by developing resilient landscapes and cities (Lafortezza et al. 2018). Urban ecological infrastructure and ecosystem services can increase the resilience of cities by enhancing their ability to cope with disturbance and climate change and adapt to climate and other global change. Green infrastructure and ecosystem services can be referred to as a form of insurance value. Enhancing green space and green infrastructure ensures a resilient supply of ESs. The results from this research could provide insights to urban planners and land managers for a better distribution of ESs.

Model performance and landscape configuration

Landscape configuration metrics appear to be strong determinants of ES trade-offs and synergies. At the patch level, although the model error rates ranged between 50 and 71%, the mean rate among the models (63.0%) was comparable to results found in other studies applying BBNs to environmental systems (Aalders 2008; Grafius et al. 2019). At the class level, the principal components were found to be most sensitive to MESH. The findings suggest that landscape configuration could drive ES interactions. Area, core area and aggregation metrics appeared to be the most influential metrics in determining ES trade-offs and synergies.

Greenspaces with a large patch size and a short inter-site Euclidean distance have greater potential to support species, provide ESs and have a better ecological quality (Kim and Pauleit 2005; Tian et al. 2014). CORE, CONTIG and ENN could be used to assess the influence of configuration on nutrient retention and carbon storage trade-offs at the patch level. In this study, the differences in sensitivity analysis results between ES trade-offs and synergies models could be due to the degree and pattern of fragmentation and the ES in question (Mitchell et al. 2015a, b).

At the class level the ES interactions showed strong relationships with the configuration metrics suggesting that mesh size is a key driver of urban ES interactions. The high levels in sensitivity between ES trade-offs and synergies models suggest that landscape configuration factors affect ES interactions and could be appropriate to assess the influence of landscape structure on ES trade-offs and synergies and for planning.

Implications

The findings show that landscape configuration affects the provision of ESs and could drive ES trade-offs and synergies. The modelling approach has potential value as a method of ES trade-offs and synergies prediction in complex landscapes and could be used to assess how ES trade-offs respond to variation of urban habitat configuration and thereby relevant to planning considerations and landscape-scale research. The findings could have implications in landscape planning and management. The indicators could be used for structural assessment of nutrient retention and carbon storage, habitat quality and pollinator abundance trade-offs and soil erosion and water supply synergies and in planning in similar urban areas in the UK.

Limitations

LM have been found to be sensitive to extent, land cover classes and spatial resolutions of remotely sensed data (Verhagen et al. 2016; Inkoom et al. 2018). Therefore, users of our suggested metrics must proceed with caution, since some of the metrics may be sensitive to changing spatial resolution (Inkoom et al. 2018). Further studies on spatial–temporal remote sensing data can reveal the potential of the use of LM and its contribution to the influence of landscape structure on ES provision. In addition, studies on the influence of configuration on ES bundles could provide a better understanding of the relationships between landscape configuration and ES provision and subsequently account for landscape configuration in landscape management.

Conclusion

We have demonstrated the feasibility of using a BBN modelling approach for predicting ES trade-offs and synergies. Our approach provides useful information on the sensitivity of ES trade-offs and synergies to habitat configuration features providing ecological understanding that is relevant to planning decisions and assessment of urban development impacts. The results show that landscape configuration affects ecosystem service interactions and demonstrates that landscape configuration metrics core area and mesh size are influential determinants of ecosystem service trade-offs and synergies. The findings suggest that the core set of metrics MESH, ENN_AM, CORE, ENN and CONTIG could be used as indicators to assess ES interactions in similar urban areas. Cities will need to manage a resilient supply of ESs for an enduring supply of services in urban systems affected by global environmental change.

Further research is needed to understand how climate interacts with and drives changes in urban ecosystems, and how these changes will affect the supply of ESs. Further research is needed to assess the influence of potential drivers on ES provision and green infrastructure which may provide insights into what management or policies could improve the provision of multiple ESs. Cross-scale comparative studies of the influence of landscape structure on ESs are needed in urban areas. The findings may contribute to planning and urban design and improved ecosystem management.