Virtual landscape creation
Our in silico experiment consisted in a factorial combination of landscape composition and fragmentation across three topographies. The virtual landscapes were created in Python using the NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), scikit-(van der Walt et al. 2014), NetworkX (Hagberg et al. 2008), NLMpy (Etherington et al. 2015), GDAL (GDAL/OGR contributors 2020), and RichDEM (Barnes 2016a) packages.
Topography
Each landscape had dimensions of 2000 × 2000 cells with 25 m grain and therefore 50 × 50 km extent. The topography of each landscape was initialised with a randomly directed slope with an elevation range of 0-200 m. Further topographic complexity was then added using the Perlin noise approach (Musgrave et al. 1989) to parameterise three types of topography: plains, hills, and mountains (Supplementary Table 1) each of which was replicated 10 times. Topographic depressions were then filled (Barnes 2016b) to calculate slope (Horn 1981) and to identify river channels based on D8 drainage (O'Callaghan and Mark 1984) that drained at least 3.75 km2.
Landcover
To produce landcover distributions each landscape was first divided into patches by randomly selecting 5000 cells and then creating patches as discrete Voronoi polygons around the selected cells. The mean elevation and slope of each patch was calculated to determine the suitability of each patch for six landcover types (Table 1) based on New Zealand patterns from an analysis of the national land cover data base LCBD5 ((MWLR) 2020). Proceeding in turn from the most topographically limited to least topographically limited landcover type, several suitable patches were randomly selected to act as seed patches. The distribution of landcover was then grown from these seed patches using an Eden growth process (Eden 1961) that randomly selected a suitable patch that neighboured an existing landcover patch. If no neighbouring patches were suitable, a new seed patch would be randomly selected. This landcover growth process continued until a landcover had reached its specified maximum proportion of the landscape. Any patches that were not attributed to a landcover type were classified as bare ground.
To examine the effects of changing landcover proportions, four different landcover scenarios were specified (Table 1). According to the hypothesised optimal land use composition template (Arroyo-Rodríguez et al. 2020; Smith et al. 2013), we set baseline land use at 30% intensive (intensive grassland and crops-horticulture), 40% extensive (grasslands and exotic forestry) and 30% indigenous (shrubland and forest), referred to as ‘intermediate – even’ henceforth. We then deviated from this baseline by modifying the proportion of land under intensive use to either 10% (extensive) or 40% (intermediate - intensive), keeping the relative contributions of intensive grassland and crops-horticulture constant. Consistent with recent intensification in New Zealand and other livestock farming regions, we converted extensive grassland to intensive land use while keeping semi-natural cover at 30% (20% native forest, 10% shrubland) and total forest cover at 30% (20% native, 10% exotic). The last scenario pictured upper-end intensification with 60% intensive use, reduction of semi-natural vegetation to a total of 10% and total forest to 15% (intensive henceforth). These land cover scenarios were within actual variability in crop, intensive grassland and native forest cover in New Zealand catchments.
To examine the effects of landcover fragmentation, each landcover was grown from 5, 10, 50, 100, 500, and 1000 seed patches – increasing the number of seed patches increased the fragmentation of the resulting landcover distributions. A factorial combination of six levels of seed patches and four different landcover proportions produced 24 different landcovers for each of the 10 replications of the three topography types, a total of 240 simulated landscapes.
Table 1 - Parameterisation six landcover types. Each landscape had topographic limit for suitability set as a function of both elevation and slope. The maximum landscape proportion for each landcover type was varied for four different scenarios.
Landcover
|
Patch maximum elevation
|
Patch maximum slope
|
Maximum proportion of landscape
|
|
|
|
Extensive
|
Intermediate
- even
|
Intermediate
- intensive
|
Intensive
|
Crops-horticulture
|
500
|
10
|
0.001
|
0.01
|
0.04
|
0.06
|
Intensive grass
|
1000
|
20
|
0.099
|
0.29
|
0.36
|
0.54
|
Extensive grass
|
2000
|
90
|
0.5
|
0.3
|
0.2
|
0.2
|
Shrubland
|
1750
|
90
|
0.1
|
0.1
|
0.1
|
0.05
|
Exotic forest
|
1250
|
50
|
0.1
|
0.1
|
0.1
|
0.1
|
Native forest
|
1250
|
90
|
0.2
|
0.2
|
0.2
|
0.05
|
Ecosystem service models
We modelled six ecosystem services essential to climate adaptation and mitigation: pollination, erosion and water quality regulation for sustainable and climate-resilient (e.g. changed rainfall) food production, carbon sequestration and greenhouse gas mitigation and landscape appreciation for recreation as a critical component of rural quality of life. Spatially insensitive (e.g. carbon stocks), proximity-based (e.g. pollination, landscape aesthetic value) and flow-based (e.g. regulation of erosion or water quality) ecological processes underpin ES supply. These combine with functions of social access and value delivery (Lavorel et al. 2020): non-spatial (e.g. climate mitigation), proximity- (e.g. pollination, recreation) or topographical flow-dependence (e.g. water quality regulation), and size thresholds (e.g. carbon trading schemes - though not considered here). We modelled supply capacity for ES representing this range of sensitivities, parameterised with New Zealand data (Table 2).
Table 2 – Ecosystem service models used and parameterisation. Constitutional sensitivities of ES models to different spatial components are indicated as dots. Colouring of model parameters for the 6 land cover types indicate low (red tones) to high (green tones) effects of respective ES supply potential.
Carbon stocks: The carbon stocks model is a look-up table relating the land cover type of each pixel to an estimated value of carbon stocks (Case and Ryan 2020; Mason et al. 2012; Thomas et al. 2021). The standardised landscape sum of carbon stocks across all pixels in each landscape was calculated and used for comparisons of carbon stocks between landscape fragmentation levels and proportions sets.
Greenhouse gas emissions: Likewise, the greenhouse gas emission model is a look-up table relating the land cover value of each pixel in the landscape to an estimated annual value of greenhouse gas emissions. Values were based on crop/land-cover and animal type-based emission values from Thomas et al. (2021) and were further simplified (e.g. assuming uniform stocking-rates) and adapted, to reflect a range of different land-uses. The landscape sum of greenhouse gas emissions for all pixels in each landscape was calculated and used for comparisons of greenhouse gas emissions between landscape fragmentation levels and proportions sets. For analyses of multifunctionality and pairwise interactions, we considered the ES of avoided emissions, calculated as 1 minus the standardised emission value for each pixel.
Erosion: We implemented a derivation of the USLE model for New Zealand, NZUSLE (Dymond 2010). The NZUSLE model is the product of a precipitation factor (P), a slope gradient factor Z, a slope length factor L, a soil factor K, and a vegetation factor U with an equation of the following form where Es is the mean annual erosion rate due to surficial processes (in t km− 2 yr− 1).
Es (x,y) = αP2(x,y) Z(x,y) L(x,y) K(x,y) U(x,y)
α is a constant calibrated with published surficial erosion rates (1.2 × 10− 3). The precipitation factor (P) requires an estimate of mean annual rainfall per pixel. This was estimated for our virtual landscapes with a sea surface level mean annual rainfall of 800mm for New Zealand. and a + 1mm lapse per 1m increase in altitude.
The slope gradient factor (Z) is calculated as; 0.065 + 4.56(dz/dx) + 65.41(dz/dx)2 where dz/dx is the slope gradient.
The slope length factor (L) is calculated as (λ/22)0.5 where λ is slope length in metres.
The slope length λ was calculated from the landscape DEM as presented by (Barriuso Mediavilla et al. 2017) and (Bolton et al. 1995).
We set the soil erodibility factor (K) to 0.2, a value corresponding to medium level erodibilty clay soil types in Dymond 2010 as soil types were not simulated for this study.
Finally, values for the vegetation factor (U) were parameterised for the land cover types based on expert assessment of the relative ability of each land cover to retain surface soil particles.
The landscape sum of erosion for all pixels in each landscape was calculated and used for comparisons between landscape fragmentation levels and proportions sets. For analyses of multifunctionality and pairwise interactions, we considered the ES of avoided erosion, calculated as 1 minus the standardised erosion value for each pixel.
Nitrogen retention was modelled using the nutrient delivery ratio model (NDR) of the InVEST ecosystem service modelling package (Sharp et al. 2020). This model uses a simple mass balance approach, describing the movement of a mass of nutrient through space and represents the long-term, steady-state flow of nutrients through empirical relationships. Sources of nutrient across the landscape, also called nutrient loads, are determined based on land covers and associated loading rates. Nutrients are then transported via surface flow (we chose to not model sub-surface flows). Delivery factors are computed for each pixel based on the properties of pixels belonging to the same flow path (in particular their slope and retention efficiency of the land cover). At the watershed level, the nutrient export is computed as the sum of the pixel-level contributions. The parameterised nitrogen loads and nitrogen retention efficiencies were determined by expert assessment of available data for each land cover (Davis 2014; Elliott et al. 2005; Ledgard 2014; Pärn et al. 2012; Sharp et al. 2020). Total nitrogen retention for each landscape was calculated as the proportion between total landscape nitrogen load and total landscape nitrogen export. For analyses of multifunctionality and pairwise interactions, we calculated the standardised proportion of landscape nitrogen retained.
Crop pollination: The pollination model defines land covers on the basis of their capacity to provide habitat for pollinators, and on their requirement for pollination (Maes et al. 2012; Schulp et al. 2014b). To simulate pollinator movement, a 500 m buffer is implemented around patches of land covers providing pollinator habitat, and the overlap between those areas requiring pollination and pollinator availability is calculated. Considering provisioning services of crop (e.g. rape, clover seeds etc.), horticultural (fruit, vegetable and horticultural seed production and honey production, we defined crops / horticulture and shrublands as land covers requiring pollination. For pollinator habitat, shrublands are assigned a pixel value of 1 (high quality pollinator habitat), extensive grasslands, exotic forests and native forests a value of 0.5 (medium quality pollinator habitat) and crops / horticulture and intensive grasslands a score of 0 (no pollinator habitat). After the extension of the buffer area around each land cover patch (buffer pixels taking the value of the origin patch), the degree of overlap between areas requiring pollination and areas of high and medium quality or no pollinator habitat are identified. As shrublands are both pollinator requiring and high quality pollinator habitat, these are always pollinated, and thus we considered the proportion of crops / horticulture pixels that are overlapped by the differing degree of pollinator availability that becomes the variable of interest. We calculated the % of crops without pollinators for comparisos of pollination availability between landscape fragmentation levels and proportions sets. For analyses of multifunctionality and pairwise interactions, we considered the ES of the proportion of crops covered by pollinator habitat (either high or medium quality), calculated as 1 minus the standardised proportion of crops without pollinators for each landscape.
Recreation: We focused on landscape attractiveness from the Recreation Opportunity Spectrum approach (Byczek et al. 2018). The recreation model comprises of 3 components of landscape attractiveness that are summed to provide a compound final raster of overall landscape attractiveness. The first component relevant to fishing, swimming, scenic value defines areas close to a watercourse (determined from the landscape DEM) as more attractive for recreation than areas distant from a watercourse. Buffers of 500-meter size are created around aquatic features in the landscape (rivers) and are classified as attractive for recreation (pixel value = 1), while the remainder of the landscape is classified 0. The second component classifies areas located in elevated areas of the landscape (hills, ridges) as being more attractive for recreation than non-elevated areas. Such areas are considered to provide scenic views attractive for recreation. A topographical position index is calculated for each DEM pixel using a circle of 200 meters diameter, and a breakpoint of 1.5 used to classify pixels as elevated (> 1.5) or not (< 1.5). Elevated pixels are given a score of 1, and non-elevated pixels a score of 0. The third component weights each pixel depending on the attractiveness of its land cover and the size of patch of that land cover. Land covers such as native forests are considered as more attractive for recreation than, for example, intensive grasslands (Richards and Lavorel). Larger contiguous areas of any given land cover are considered as more attractive for recreation than small, disjunct areas of that land cover (Cordingley et al. 2015). Each pixel is assigned a subjective recreation attractiveness score according to its land cover (Table 1). The model then calculates the proportion of the landscape occupied by each patch and multiplies the land cover attractiveness score of each pixel by the proportional size of the patch in which the pixel is found. Finally, the three components of landscape recreation attractiveness are summed to provide a compound score of recreation attractiveness. The mean of this value across each landscape is used for comparisons of recreation attractiveness between landscape fragmentation levels and proportions sets. For analyses of multifunctionality and pairwise interactions, we calculated the standardised value of the recreation attractiveness score.
The six ES models were coded into a Snakemake pipeline (Mölder et al. 2021), using: Python nad Python packages RasterStats (Perry 2021), Pandas (The pandas development team, 2021), SciPy (Virtanen et al. 2020), InVEST (Sharp et al. 2020), Rasterio (Gillies et al. 2013); R and R packages landscapemetrics (Hesselbarth et al. 2019), raster (Hijmans 2021), and dplyr (Wickham et al. 2021); GDAL (GDAL/OGR contributors 2021); GRASS GIS (Neteler 2021); Docker and Debian Buster.
Data analysis
Landscape pattern was quantified using common spatial metrics. After screening for correlations we retained: mean and land cover specific patch number and mean patch size, nearest neighbour distance, connectivity and Shannon diversity for characterising individual landscapes (Rieb and Bennett 2020).
Responses of individual ES to combinations of landscape composition and fragmentation were tested using General Linear Models with topography (3 types), and the combination of intensification level (categorical land cover parameter sets) and mean patch size (continuous) as explanatory variables.
We quantified multifunctionality across the six ES, and again for only pollination, nitrogen retention and recreation, using four indicators. The first two indicators consider the magnitude and evenness of provision by each virtual landscape as a whole, while the second two consider within-landscape variation by comparing scores across different landscape units. The four indicators were:
-
Multifunctionality at the landscape scale, quantified as the mean sum of scaled provision for 6 (resp. 3) ES across the landscape
-
Evenness of provision at the landscape scale, quantified as the Pielou evenness (J) of scaled provision scores for the 6 (resp. 3) ES across the landscape. L ow values indicate specialisation towards particular ES, while high values indicate even (though not necessarily high) performance across ES
-
Mean landscape scale multifunctionality was quantified as the mean multifunctionality scores for all pixels within a virtual landscape. Multifunctionality within each pixel was quantified as the mean sum of scaled provision for the 6 (resp. 3) ES of individual pixels
-
Heterogeneity in landscape scale multifunctionality was quantified as the standard deviation in multifunctionality scores for all pixels within a virtual landscape. Multifunctionality within each pixel was quantified as the mean sum of scaled provision for the 6 (resp. 3) ES of individual pixels
For the landscape scale indicators, each ES score for each landscape was scaled against the minimum and maximum values across all landscapes. For the patch scale indicators, each ES score was scaled against the minimum and maximum values found in any one pixel across all virtual landscapes.
We modelled all multifunctionality indicators as mixed-effects generalized linear models fitted using penalized quasi-likelihood and assuming a quasibinomial error structure. Models were fitted using the topography type factor as a random effect and the mean patch size and landscape composition treatment group as fixed effects. An interaction was fitted between the two fixed effects. Models were fitted using the glmmPQL R function (Venables and Ripley 2002).
Trade-offs and synergies between pairs of ecosystem services were quantified following the method proposed by (Vallet et al. 2018). This approach was previously used to analyse a cloud of points resulting from a series of iterative but not optimised scenarios – similar to our case (Vallet et al. 2018). For each pair of ecosystem services, the scaled indicator scores were plotted against each other, and the overall shape of the bounding box enclosing the cloud of points was characterised as the alpha-convex hull (Pateiro-Lopez and Rodriguez-Casal 2019). The shape of the bounding box was quantified using the shape index I, which characterises the “roundness” or “elongatedness” of the cloud (Vallet et al. 2018). The correlation coefficient r of the ES scores was quantified as an indicator of the overall trend in the cloud of points (Vallet et al. 2018). The Pareto frontier of each pair of ecosystem services was quantified using non-dominated sorting of the scaled scores.
The shape of the cloud of points representing each trade-off or synergy relationship may be influenced by the input parameters used to generate the virtual landscapes. To explore the contribution of the fragmentation and landscape composition treatments in deriving the shape of the trade-off or synergy, we conducted a resampling sensitivity analysis. For both the landscape composition and fragmentation treatments, we conducted a series of simulations in which a randomly selected treatment group was removed from the dataset, and the remaining datapoints were resampled with replacement to make up a dataset of equivalent size to the original (n = 720). For each replicate we conducted all pairwise ecosystem service comparisons, and the resulting clouds of points were analysed to quantify the shape parameter I and the correlation coefficient r of the overall cloud. The resampling was conducted 1000 times, and 95 percentile ranges of the parameters were calculated to quantify the sensitivity to fragmentation and landscape composition. All statistical analyses were conducted in R version 4.0.4. (R Core Team, 2021).