Copyright © 2016 by the author(s). Published here under license by the Resilience Alliance. Meacham, M., C. Queiroz, A. V. Norström, and G. D. Peterson. 2016. Social-ecological drivers of multiple ecosystem services: what variables explain patterns of ecosystem services across the Norrström drainage basin?. Ecology and Society 21(1): 14. http://dx.doi. org/10.5751/ES-08077-210114 E&S Research, part of a Special Feature on Programme on Ecosystem Change and Society (PECS): Knowledge for Sustainable Stewardship of Social-ecological Systems Social-ecological drivers of multiple ecosystem services: what variables explain patterns of ecosystem services across the Norrstrom drainage basin? Megan Meacham'. Cibele Queiroz1'2. Albert V. Norstrom' and Garry D. Peterson' ABSTRACT. In human dominated landscapes many diverse, and often antagonistic, human activities are intentionally and inadvertently determining the supply of various ecosystem services. Understanding how different social and ecological factors shape the availability of ecosystem services is essential for fair and effective policy and management. In this paper, we evaluate how well alternative social-ecological models of human impact on ecosystems explain patterns of 16 ecosystem services (ES) across the 62 municipalities of the Norrstrom drainage basin in Sweden. We test four models of human impact on ecosystems, land use, ecological modernization, ecological footprint, and location theory, and test their ability to predict both individual ES and bundles of ES. We find that different models do best to predict different types of individual ES. Land use is the best model for predicting provisioning services, standing water quality, biodiversity appreciation, and cross-country skiing, while other models work better for the remaining services. However, this range of models is not able to predict some of the cultural ES. ES bundles are predicted worse than individual ES by these models, but provide a clear picture of variation in multiple ecosystem services based on limited information. Based on our results, we offer suggestions on how social-ecological modeling and assessments of ecosystems can be further developed. KeyWords: ecological footprint; ecologicalmodernization; ecosystem service bundles; land use change; location theory; Stockholm; Sweden INTRODUCTION Ecosystem services are coproduced by humans and nature as a result of interactions between ecological functions and societal management and demand (Schröter et al. 2012, Reyers et al. 2013). However, understanding how this social-ecological coproduction of multiple services occurs in specific landscapes is poorly understood. Social-ecological interactions often produce distinct patterns of ecosystem services across landscapes, in the form of coherent sets of ecosystem service bundles (Bennett et al. 2009, Raudsepp-Hearne et al. 2010, Hanspach et al. 2014, Queiroz et al. 2015). Management or policy measures that aim to enhance a single, particular service will miss this complexity and can lead to perverse effects caused by hidden trade-offs among services such as those between crop yield and water quality (Tilman et al. 2002, Zhang et al. 2007), and timber extraction and carbon storage in forests (Putz and Romero 2001, Nelson et al. 2008). Trade-offs are often found between provisioning servicers and regulating or cultural services. Explicitly focusing on ecosystem service interactions and, in particular, bundles of ecosystem services is useful for the management of complex landscapes because it avoids those pitfalls. Furthermore, it helps identifying interventions that can have simultaneously desired effects on multiple ecosystem services (Queiroz et al. 2015). Thus, a crucial challenge is to identify and understand the key drivers, and their interactions, that produce these coherent sets of ecosystem service bundles across landscapes. Because governments, NGOs, and companies are increasingly interested in incorporating ecosystem services into their decision making, there is a need for ways of predicting and quantifying the production, management, and use of ecosystem services. However, the current available tools are often complex and do not match data availability or existing decision or policy scales (Daily et al. 2009, de Groot et al. 2010). We test how easily accessible data can be used to understand the distribution of multiple ecosystem services across a spatial context that is relevant for decision making, e.g., multiple municipalities in a human-dominated landscape. Previous research has shown that spatial patterns of provisioning, regulating, and cultural ecosystem services produced in a landscape are driven by a mix of social factors, such as distance from a city, and ecological factors, such as topography and land use (Raudsepp-Hearne et al. 2010). We compare the ability of different models of human-environment interactions to explain patterns of ecosystem services observed across the Norrstrom drainage basin. Different academic disciplines have emphasized distinct variables in their analysis of human-environment interactions (Briassoulis 2000, Lambin et al. 2006, Hersperger et al. 2010). We use established theories of human environmental interaction from human ecology, political science and economics, ecology, and geography (Table 1). We adapt four dominant alternative models of human-environment interactions to assess the importance of different potential drivers of ecosystem services as well as understand the strengths and limitations of these theories in the context of ecosystem services. We define drivers as the fast and slow changing variables, including exogenous controls, like slope (Carpenter et al. 2009). Ecologists have proposed that land use or land cover integrates social-ecological interactions and is therefore a useful proxy for assessing how humans interact with landscapes. This proposition is less of a model than an assumption that has been underlying many ecosystem services assessments, in which land use or land cover is accepted as the main determinant for the potential supply of ecosystem services (Costanza et al. 1997, Burkhard et al. 2009, Nelson et al. 2009). Political science and economics have proposed ecological modernization models that suggest that environmental quality improves as societies become wealthier and more developed (Kuznets 1955, Grossman and Krueger 1995, Mol 2003). From human ecology, the ecological footprint approach has highlighted the importance of social driver, e.g., population density, affluence, and technology, in driving negative environmental impact and 'Stockholm Resilience Centre, Stockholm University, 2The Beijer Institute of Ecological Economics, Swedish Royal Academy of Sciences Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Table 1. Models used to estimate ecosystem services and the specific drivers and proxies used in each model. Model Drivers Proxy measure Units Land use Arable land Proportion of arable land % Grazing land Proportion of grazing land % Forest area Proportion of forest % Built up land Proportion of built up land % Ecological modernization Education Adults with some secondary education % Wealth Difference between total assets and total debts per capita SEK Ecological footprint Income Net annual income SEK/year Population density Log Population density People/km2 Location theory Distance from Stockholm Distance from Stockholm Meters Slope Average slope Meters consequently affecting ecosystem service distribution patterns (Stern et al. 1992, Chertow 2000, Harrison and Pearce 2000, York et al. 2003). Finally, geography has shown how distance and the cost of transport strongly shapes land use (von Thiinen 1826, Alonso 1964, Griffith 1999), suggesting that variables such as distance from a city and topography, have a strong impact on activities performed across the landscape. METHODS Study area and ecosystem services Our study region is the Norrstrom drainage basin, which is situated in southcentral Sweden. It covers 22,650 km2 and includes two of Sweden's largest lakes, Lake Hjalmaren and Lake Malaren as well as Stockholm, Sweden's largest city. The Stockholm metropolitan area, which is inhabited by approximately 1.5 million people, occupies the eastern part of the region and contains most of the population and economic activity of the region, and is the economic and political center of Sweden. Lake Malaren supplies the drinking water for the Stockholm metropolitan area. The Stockholm region has a maritime climate; it is urbanized, and has high population growth. Around Lake Malaren, the region is dominated by agricultural land, which is primarily used for the cultivation of cereals and rapeseed, and pasturage. The westernmost parts of the region have a more continental climate. Commercial forests characterize the land cover and land use in this part of the basin. This area has been experiencing a decline in population over the last 30 years as people have moved to more urban areas. A variety of ecosystem services are produced, managed, and used in the Norrstrom drainage basin. Not only do the lakes in this region provide drinking water for the surrounding population, but they also serve as an area for recreation as well as a reservoir, retaining nutrients before they are released into the Baltic Sea. Pollinator habitats are abundant around agricultural areas, supporting insect pollinated crops and pasture areas. In Sweden, there is demand for summer cottages, outdoor recreation, the reporting of sighted biodiversity, and moose hunting, all these supported by both open habitats and forests. The forest areas across the basin also provide a range of commercial products, such as timber, paper fiber, wild berries, and mushrooms among others. In a previous study we quantified and mapped 16 ecosystem services (Queiroz et al. 2015), which we use as the ecosystem service information in this study (see Appendix 1). The services included were those that had been expressed as of interest by policy documents, regional planners, and other stakeholders in the region. Table 2 provides a complete list of those ecosystem services. Queiroz et al. (2015) showed that the Norrstrom drainage basin is characterized by five ecosystem service bundles, which are spatially clustered across the basin. Each bundle type is characterized by a distinctive arrangement of the ecosystem services, represented by flower diagrams in Figure 1. Two of these bundles, Mosaic cropland horse and Mosaic cropland livestock, are characterized by high values of regulating services and provisioning services such as crop production, but differ in their levels of livestock production. The Forest and towns and Remote forests bundles are defined by the high provision of forest products, but differ with respect to certain cultural services such as moose hunting. The final bundle, Urban, is spatially constrained to the municipalities in metropolitan Stockholm and had the highest values of cultural services such as cross-country skiing and biodiversity appreciation. Selecting drivers and model creation Four candidate models, Land Use, Ecological Modernization, Ecological Footprint, and Location Theory, were created by carefully considering different key theories of human-environmental interactions and distilling the different driver variables that each theory emphasizes. The theories included in this paper are studied and cited theories on human-environment interactions prominent in different academic disciplines (Dietz et al. 2007). Each theory postulates what are key driving forces behind these interactions. We interpreted these driving forces into measurable driver variables to construct each model. We tested how these four candidate models predicted patterns of ecosystem service production in the study region. We also created a reference model that combines all the drivers used in the different candidate models. We used information from public databases to quantify the different social, ecological, and geographic variables that the alternative key theories suggested were drivers of human nature interactions (see Appendix 2). Variables used in the candidate models (Table 1), to capture these drivers, were also based on expert knowledge of the region, data availability, and spatial heterogeneity across the landscape. Below we describe the candidate models, and explain the driver variables included, in more detail. Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Table 2. Ecosystem services by category and the units in which they were measured. Ecosystem service Unit Provisioning Wheat production km2 Cattle Cattle/km2 Pig production Pigs/km2 Sheep Sheep/km Forest products km2 Regulating Crop pollination km2 Nitrogen retention (N retention) Average fraction of l-(net nutrient load/Gross nutrient load) Phosphorous retention (P retention) Average fraction of l-(net nutrient load/Gross nutrient load) Standing water quality Average water quality/km2 Running water quality Average water quality/km2 Cultural Moose hunting Kills/km2 Outdoor recreation km2 Summer cottages km2 Horse back riding Horses/km2 Cross country skiing Stations/km2 Biodiversity appreciation Species reported/km2 Fig. 1. Bundles of ecosystem services in the Norrstrom drainage basin (adapted from Queiroz et al. 2015). Each ecosystem service bundle is represented by a flower diagram, in which the area of each wedge represents the quantity of an ecosystem service. The names of the bundles are descriptive of the social-ecological context and all municipalities belonging to a particular bundle are shaded. Land use (arable land, grazing land, forest area, built-up land) The land-use driver category represents different suites of ecological factors that are managed for. The use of land use stems from human ecology and emphasizes that it is the biophysical landscape that constrains human activity (York et al. 2003). Land use has been widely used as a proxy of ecosystem service supply in ecosystem service assessments (Costanza et al. 1997, Burkhard et al. 2009, Nelson et al. 2009). Ecological modernization model (wealth, education) Ecological modernization predicts that increasing socioeconomic development results in ecological degradation until a point when environmental conditions improve as societies become increasingly mature and affluent and begin to demand environmental quality (Kuznets 1955, Grossman and Krueger 1995). The ecological modernization perspective is sometimes represented by the environmental Kuznets curve and postulates that economic development will provide solutions to environmental problems (York et al. 2003). Ecological footprint (population density, income) The ecological footprint concept emphasizes that the environment provides resources, living space, and acts as a sink for waste. As population density and economic intensity distances populations from the supply of resources and their waste, there is inefficient use, degradation, and exploitation of these regions (Jorgenson 2003, York et al. 2003). The ecological footprint concept is an articulation of the IPAT formulation (Stern et al. 1992, Chertow 2000, Harrison and Pearce 2000) that states that environmental impact is equal to the multiplication of population density, affluence, and technology. Location theory model (distance from Stockholm, slope) Location theory stems from geography and economics and supposes that the distance from an urban center will determine the land-use type and activities available for that area. The values of land, along with the heterogeneity of the landscape are the constraints in the spatial distribution of activities in the landscape (von Thiinen 1826, Alonso 1964, Griffith 1999). Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ We calculated values for each of these drivers for the 62 municipalities in this landscape. Some of the drivers are correlated with one another, for example distance from Stockholm and income, but these correlated drivers were included in the analysis because of their different theoretical contributions, for example the Location Theory and Ecological Modernization models. Table 1 lists the drivers used across the four models tested: Land Use, Ecological Modernization, Ecological Footprint, and Location Theory. Data analysis Mapping of driver variables We used QGIS 2.0.1 to map the spatial distribution of each individual driver by municipality and R to conduct all statistical analyses and produce all the figures (R Core Team 2013). The R packages ggplot2 (Wickham 2009), sp, and maptools (Bivand et al. 2013) were used to create the figures. Relationship between driver variables We identified correlations among drivers by comparing all the drivers to each other using the Pearson's correlation coefficient. We compared how the alternative candidate models predicted individual ecosystem services by fitting ordinary least-squares multiple linear regression models of each model to each service, and then compared these modes using adjusted R2 and the Akaike information criterion (AIC) to compare models for predictive power and parsimony. When comparing the fit of several models, AIC provides a criterion (penalized log likelihood) with which to find the model that best explains the data with a minimum of free parameters (Akaike 1974, Burnham and Anderson 2004). We also used a correction for small sample sizes (AICc). Relationship between drivers, candidate models, and ecosystem services We conducted a principal component analysis of instrumental variables (PCAIV), using the R package ade4 (Dray and Dufour 2007), to determine how the assessed drivers and assessed ecosystem services were related to one another. PCAIV is a method, similar to redundancy analysis and canonical correspondence analysis, that is used in ecology to study the relationships between the composition of species communities and their environment, by matching a sites-by-environmental-variables table and a sites-by-species table. Here we adapted it to analyze the relationship between sites-by-drivers and sites-by-ecosystem services. To identify the relative importance of each individual driver variable (independent of correlation among variables) in predicting individual ecosystem services, we used hierarchical partitioning. This is a statistical method that analyzes all possible models in a multiple regression to identify the relative contribution of each variable to the total variance, both independently and in conjunction with the other variables, to estimate the impact of each variable (Chevan and Sutherland 1991, MacNally 2002). Specifically, we used the hierarchical partitioning to test the independent effect of all drivers, except slope on the prediction of each ecosystem service. We excluded slope, because hierarchical partitioning is only able to simultaneously test 9 variables, and slope was the only variable that was not significant in any model of ecosystem services. We used the R package "hier.part" (MacNally and Walsh 2004). We assumed Gaussian errors and calculated goodness of fit using R2, and the statistical significance of the independent effect of each variable was estimated using a randomization approach (n = 1000) to calculate Z scores (MacNally 2002). We used random forests, using the R package randomForest (Liaw and Wiener 2002), to assess and compare how well the alternative candidate models predicted both individual services and bundles of ecosystem services across the study landscape. This approach provides robust results for highly correlated predictor variables and small sample sizes (Prasad et al. 2006, Cutler et al. 2007). The random forest routine was used to predict the most likely ecosystem service bundle in each of the municipalities, but also provided the votes the random forest gives for each bundle. This allows the uncertainty of different predictions to be compared. The random forest analysis constructs a multitude of decision trees and produces a mean prediction of the individual trees. By comparing across individual trees the random forest allows the relative importance of different predictor variables to be assessed. We did this by calculating 100,000 trees for each model, and then calculating the accuracy of the bundle prediction as the percentage of correct predictions; because some ecosystem service bundles are more similar to each other we also assessed how well the prediction of bundles actually predicted ecosystem services in each municipality. We also used the results of the random forest to predict ecosystem services in each bundle by two methods. First, the average level of the set of ecosystem services predicted determines the bundle for each municipality. The second method is the voting method. Each iteration of the random forest selects a bundle and that is counted as one vote for that bundle and the weighted average of the votes for different bundles of ecosystem services is the result. These two methods only produce different results if there is significant vote for different bundles of ecosystem services. Comparing the results of the linear models to the results of the randomForest is not straightforward, however, we can compare the predictions of ecosystem services from each approach by calculating R2. This allows us to compare model predictions to observations, but this approach does not adjust for the increased complexity of the linear models, and the risk that they overfit the data. More robustly comparing model predicts require larger numbers of municipalities or other ecosystem service delivery units, and represents an area of future research. To test for unaccounted effects of spatial autocorrelation on our analysis we calculated Moran's I on the residuals of the model's predictions of ecosystem services, accounting for shared boundaries between municipalities. We also plotted model residuals for each municipality by ecosystem service and model. This approach allows us to assess whether models have adequately accounted for spatial autocorrelation among the ecosystem services, or if there is still substantial unexplained autocorrelation, and we use it as diagnostic of how well models explain ecosystem service pattern rather than to create spatial regression models (Dormann et al. 2007). RESULTS Distribution of individual drivers Each of the 10 drivers analyzed showed distinct patterns of spatial distribution across the basin (Fig. 2). All social drivers, wealth, income, population density, and education, exhibited highest Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Fig. 2. Distribution of drivers across the Norrstrom drainage basin's municipalities. The values of drivers are shown in equal interval quintiles. Darker shades of blue represent higher values for the drivers. Levels of the drivers are shown in the legends for each driver. Population Density Wealth Income Education Arable Land Leg Population Density □ 1-2.Ö □ 2.B-4.2 ■ 4-2-5lS ■ 5.3-7.4 ■ 7.4- B ■ No Data Log Wealth haarte Education □ 12 - 12.B □ 187- 1B4.B □ 11-17.2 □ 12.9 - 13.6 □ 1B4.6-2222 □ 17.2-214 ■ 13.6 - 14.4 □ 2222-249.B ■ 234-29 B ■ 14.4 - 16.2 ■ 249.B - 277.4 ■ 29.6-35.3 ■ 16.2 - 16 ■ 277.4-305 ■ 35.3-42 ■ Ha U ■ No Data Built up Land Forest Distance from Stockholm □ 0-61.4 □ 61.4-122.B ■ 122.B-1B4.2 ■ 1S4.2-245.B ■ 245.9-3D7 values in the Stockholm metropolitan area. In contrast, land-use drivers exhibited more complex patterns. For example, forest area had relatively high values across the basin with the exception of the urban municipalities around Stockholm, while the percentage of built-up areas was highest in the Stockholm municipalities. The proportion of arable and grazing land was relatively higher on the municipalities around lakes Malaren and Hjalmaren and some coastal municipalities north of Stockholm. Finally, the geographical drivers highlight the contrast between the northwest part of the basin with higher slopes, located further away from the city of Stockholm, and the relatively flat areas around lakes Malaren and Hjalmaren, that were also closer to Stockholm. Interactions between individual drivers We found strong negative and positive correlations between several of the drivers. Notably, the strongest correlations were not found between drivers of the same category (social, ecological, or geographical). Instead, two distinct clusters of positive correlations emerged; the first capturing the positive relationship between forest area and distance from Stockholm, while the second showed positive correlations between wealth, income, population density, education and built-up area (see Fig. A3.1). These two groupings are negatively correlated with each other; distance from Stockholm and forest area are both negatively correlated with wealth, income, population density, education and built-up area. Arable and grazing land were strongly positively correlated with each other and exhibited a strong negative relationship with slope and a weaker negative correlation with forest area. Grazing areas also showed a strong negative correlation with distance from Stockholm, whereas arable land showed no strong correlation with this variable. This suggests a trade-off between the forested areas and urban type areas. These types of relationships are explored more closely in the PCAIV analysis that follows. Using drivers to predict the distribution of individual ecosystem services The relationship between drivers and individual ecosystem services is represented by the results obtained with the PCAIV (Fig. 3). Figure 3a represents the distribution of ecosystem services with relation to the drivers. Figures 3b and 3c show the distribution of municipalities based on the relationship of their ecosystem services and drivers. The lighter colors in Figure 3b show the municipalities where the driver variables better explain the ecosystem services. The different colors in Figure 3c represent the different bundles of ecosystem services to which the municipalities belong. This analysis identified two dominant axes of the relationship between drivers and ecosystem services that together explain about 75% of the variation across municipalities. The first axis represents urbanism or connection to Stockholm; it is defined by a combination of population density, wealth, income, education, and built-up land. Thus, it separated urban Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Fig. 3. Principal component analysis of instrumental variables (PCAIV) of ecosystem services with respect to drivers as well as variance among municipalities. The X axis corresponds to rural to urban gradient (left^right), while the Y axis is shows a gradient between agricultural and forest land (bottom^top). The three plots represent (a) Distribution of individual ES and drivers, (b) Distribution of the 62 municipalities of the Norrstrom drainage basin, (c) Distribution of the ecosystem service bundles across these municipalities. p\ Ludvíka ■ - ' , Siqluns Ljusnarsberg 3 SKJnnskattebera Hec Laxa Fa( NoraNoit Vansbro imora srsta Gaqnef erga Degerfors Sti Hallslah SollenWÄ ündesberg Askersund Hábo jjy Östhamri a nam nidi immar Tierp ^ala Borlänge Huddinge eA^s\ad'ebacken Stockholm ™mBotkyrka Solná ř Sundbyberg acka „ ^Kopinq LekebersUp[)sa|a Up Plen Knivsta EskilstunaHeb Karnneholm gäte Örebro Enkoping Up Vingaker Kunqsör NoiTtälje planÜsVäsby ^ jntuna i-allefo's aby Järfälla Danderyd Södertälje ■ ^nds-Bro H umla Ekera a D a 2.5 5 Axisl D medium • high c) Ludvika •Sjgtun; Ljusnarsbt Skinnsj^neBei Degeflbr Hallstahamjjar Karlskoc a Sollen1 ■" Lindt Asfcersur ij Vansbro led! mora TgFr ^fagnef h orberg rssiit ihammar I rierp 5fl Ö*iar Ajpuya Stock iJlkvařrí0'^, L . íarimar Bo'kyrka ýolna Gnesta Köping VasfPras Lel&berg Uppsal; Val intuna yen Kcuvsta KatrinefWm n Cfebro ^ Eloping UF^nds.Bro Vřftjáker Kujigsi Sundbyberg Nicka* Norrtälje . äs íplařids Väsbi ml lällefors Mms Dandemj SödelSe ES Bundles • Mosaic Cropland Livestock • Mosaic Cropland Horse t Remote Forest a Forest and Towns • Urban municipalities located near Stockholm, and where population density, education, wealth, and income are relatively higher, from more rural municipalities. The second axis represented a gradient from agricultural land to forest dominated land. This axis was defined primarily by land use, distance from Stockholm, and slope. It differentiated forested municipalities with higher slopes of the Northwest part of the basin from the flatter agricultural land located around lakes Malaren and Hjalmaren (Fig. 3c). Without considering covariation among ecosystem services, the drivers explain 59% of the variation among ecosystem services across municipalities. To see the random forest analysis ranking of all the drivers see Figure A4.1 There were large differences regarding how well the different candidate models predicted individual ecosystem services. (Fig. 4). Some services, e.g., forest products, cross-country skiing, and biodiversity appreciation, were strongly predicted by all candidate models, including the reference model. Others, e.g., wheat production, cattle, and pollination, were predicted well by one candidate model, Land Use, and the reference model. Finally, a number of services, e.g., pigs, N-retention, P-retention, outdoor recreation, moose hunting, and summer cottages, were poorly predicted by all five models. Fig. 4. Best models for predicting individual ecosystem services. Each bar graph represents each model's R2 value for predicting each ecosystem service. ecol footprint location theory ecol footprint location theory ecol footprint location theory Forest Products Cattle ■ i i ■ ■ ■ Wheal Pollination N Retention P Retention ■ ■ Standing Water Quality Running Water Quality Moose Hunting Outdoor Recreation ■ 5 ■ 1 1 ■■■ ■■■ Summer Cottages Horseback Riding Cross Country Ski Biodiversity Appr'n land use - | ecol modernization-ecol footprint -location theory - | 0 0.25 0.50 0.75 1 0 0.25 0.50 0.75 1 0 0.25 0.50 0.75 1 0 0.25 0.50 0.75 1 AdjR2 Overall, our results show that the reference model performed best in predicting the different ecosystem services. Land Use was the candidate model that best predicted individual ecosystem services (highest adjusted R2 for 8 out of 16 services). Specifically, this model had the highest performance for predicting all provisioning services assessed in this study, three of the regulating services (pollination, standing water quality and running water quality), and one cultural service (cross-country skiing). Nevertheless, some regulating services and the maj ority of cultural services were clearly not well predicted by land use drivers. The Location Theory model was best predicting N and P retention and horseback riding, although these were better predicted by the reference model. The two candidate models that were based on Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Fig. 5. Relative ability of the models to predict individual services and bundles within municipalities across the Norrstrom drainage basin. 5a represents the gradient of accuracy in predicting the variation of individual ecosystem services within each of the 62 municipalities. The R2 of predicted versus observed ecosystem services in each municipality is shaded from dark green (perfectly correct), to opposite (dark red). 5b shows which municipalities the models correctly (match) and incorrectly (mismatch) predicted the ecosystem service bundle. 5c represents the gradient of accuracy by each model by voting for bundle type (its weighted average of random forest support for each bundle multiplied by the average level of ecosystem services in each bundle). 5a ESmodels location theory R2= 0.64 ESmodels ecological footprint R2= 0.61 ESmodels ecological modernization R2= 0.6 ESmodels land use ESmodels all £»*fe *Sddt 5^ PPr ppr pHr fST Municipality R2 ■ -1 - -0.8 ■ -0.8 - -0.6 ■ -0.6 - -0.3 □ -0.3--0.1 □ -0.1 -0.1 □ 0.1 -0.3 □ 0.3-0.6 ■ 0.6-0.8 ■ 0.8-1 5b location theory Accuracy= 0.53 ecological footprint Accuracy= 0.58 ecological modernization Accuracy= 0.47 land use Accuracy= 0.73 all Accuracy= 0.73 ■ mismatch ■ match 5c BundleVotes location theory R2= 0.56 BundleVotes ecological footprint R2= 0.58 BundleVotes ecological modernization R2= 0.52 BundleVotes land use R2= 0.53 BundleVotes all R2= 0.55 Municipality R2 ■ -1 - -0.8 ■ -0.8--0.6 ■ -0.6--0.3 □ -0.3--0.1 □ -0.1 -0.1 □ 0.1 -0.3 □ 0.3-0.6 ■ 0.6-0.8 ■ 0.8-1 socioeconomical drivers, Ecological Modernization and Ecological Footprint, were only able to better predict a few services. The Ecological Modernization model predicted the distribution of moose hunting and summer cottages relatively better, while the Ecological Footprint model predicted biodiversity appreciation better. Importantly, while relatively better predicted by the Ecological Modernization model, summer cottages and moose hunting were not well predicted by any of the sets of drivers included in this study. The overall capacity of the different models for predicting individual ecosystem services within each municipality is shown in Figure 5, given by an R2 by municipality (darker tones of green indicate higher R2 and thus a better performance of the model). An overall R2 for all municipalities is also shown in the figure. Although the reference model had the highest ability to predict the distribution of individual services (R2 = 0.79), the Land Use model performed almost as well (R2 = 0.72). The Ecological Modernization model and the Ecological Footprint model had the lowest capacity for explaining the data variability (R2 = 0.60 and R2 = 0.61, respectively). The figure shows that none of the assessed models had the ability to predict the distribution of individual services well in the urban municipalities close to Stockholm. The hierarchical partitioning analysis showed that the ability to predict ecosystem services was shared among variables, and that different ecosystem services were best predicted by different variables (see Fig. A5.1). Land variables did well at predicting agricultural ecosystem services, and log population density, log wealth, and distance to Stockholm did well at predicting many other ecosystem services. No variables predicted outdoor recreation, and with the exception of built-up land being a good predictor of cross-country skiing, many variables were weak predictors of most cultural services. Using drivers to predict bundles of ecosystem services We used the same candidate models to predict the spatial distribution of ecosystem service bundles across the basin (Fig. 5). The results for the variation of bundles predicted by the different models, through two different methods, are shown as Figs. 5b and 5c. Although the first (Fig. 5b) only measured if bundles were rightly or wrongly predicted in each of the municipalities, the second gave a measure of how much of the bundles variation within each municipality is explained by the different models, generating a gradient of accuracy from low to high predictability (Fig. 5c). The Land Use model equaled the reference model for the highest accuracy in predicting the match or mismatch between the predicted and the actual distribution of Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ bundles across the basin (R2 = 0.73), a significantly better performance than the other three models (Fig. 5b). All tested models had a low ability for predicting the differentiation between the distribution of the two agricultural bundles, Mosaic crop-livestock and Mosaic crop-horses, which included most municipalities around lakes Malaren and Hjalmaren. These municipalities were also the ones that presented higher multifunctionality, which indicates higher complexity (Fig. 1; Queiroz et al. 2015). When considering the gradient of accuracy (from low to high predictability depending on the amount of bundles variation explained by each model), the models based on socioeconomic variables performed better than the Land Use model (Fig. 5c). The Ecological Footprint model did best (R2 = 0.58) followed by Location Theory (R2 = 0.56). These both performed better than the reference model (R2 = 0.55). Both methods showed that none of the models were able to predict the distribution of ecosystem service bundles for the Stockholm municipalities well (urban bundle group), which had a high expression of cultural services. Unexplained spatial autocorrelation was found among the residuals for almost all ecosystem services for all models, for both bundles and individual ecosystem services. The pervasiveness of autocorrelation suggests that there are spatial patterns in the region and that the models are not accounting for. These are exploratory models are they are not able to account for all the spatial variation. This indicates that other drivers are needed to be included in the models or there needs to be some type of spatial interaction included. DISCUSSION The distribution of both individual services and bundles of ecosystem services across a human-dominated landscape can be well predicted by sets of relatively few drivers based on publicly available data (Fig. 5). There is a cost to data, and the ability to create understanding about the ecosystem service context using limited variables is important particularly in data scarce regions. We found that the reference model that combined all drivers, using social, ecological, and geographical variables together, performed better than the models based on one single theory or set of drivers alone. These findings suggest that ecosystem service assessments could be improved by adopting a more social-ecological approach that integrates multiple social, geographical, and ecological theories, and utilizes diverse types of data. Several of our results suggest that integrating social-ecological interactions may simplify rather than complicate the analysis of ecosystem services. First, we found that by integrating different models and variables we can predict a substantial amount of the observed variation among services. Second, our results showed the existence of strong correlations among different types of drivers, such as the correlation between income, wealth, built-up land, and distance to Stockholm. This suggests that different complementary sets of drivers can capture large proportions of social-ecological variation. Third, using ecosystem service bundles to predict average values of individual services was enough to capture a substantial variation of ecosystem service patterns across the landscape. These three findings suggest that new social-ecological models of how different types of social, geographic, and ecological drivers are related to the distribution of ecosystem services may offer rapid and reasonably accurate ways to assess their distribution across landscapes. However, although predicting the distribution of individual ecosystem services using bundles worked reasonably well, this was not the case in urban areas. Cultural ecosystem services are one of the defining features of urban areas, and the urban bundle in this landscape. One potential explanation for the inability of the models to successfully predict the bundles for these municipalities is that the drivers included in this study performed relatively poorly in predicting cultural ecosystem services. Also, the delineation between bundles of ecosystem services in the Norrstrom drainage basin is less distinct (Fig. 3c) than that of other regions where similar analyses have been done (Raudsepp-Hearne et al. 2010, Turner et al. 2014). This lack of clear boundaries between bundles is the reason behind the discrepancy between the results obtained with the method used in Fig. 5b and those obtained using the voting method described in Fig. 5c. Although the models failed to correctly predict bundles for many of the municipalities in Fig. 5b, predictions were still quite accurate with the weighted average of the votes for different bundles of ecosystem services (Fig. 5c). Predicting bundles of ecosystem services is a fairly quick and simple strategy for understanding the interaction among services and inferring something about the social-ecological context of an area. However, this method is not likely to be useful unless bundles of ecosystem services are more clearly defined within a landscape. Still, even in those cases where bundles are not very distinct, to explore the relationship between services, for example, by correlation analysis, is always a useful approach and avoids some of the problems of single service approaches. Future work should test which approach, bundles or interactions between individual services, is more useful in different landscapes. Our results partially support the use of land use as an integrated social-ecological indicator of ecosystem services, especially provisioning services, as it has been used in previous studies (e.g., Costanza and Folke 1997, Chan et al. 2006). However, we highlight certain limits to this approach. For example, the Land Use model poorly predicted a number of regulating and cultural services. Specifically, two cultural services that could be thought of as closely tied to land use, horseback riding and moose hunting, were best predicted by the Location Theory model and the Ecological Modernization model, respectively. This suggests that the reliance on land use as an indicator of ecosystem services may miss some types of regulating and cultural services. Despite success in predicting ecosystem services, there is a lack of a clear theory for explaining the production of observed patterns of ecosystem services. No individual candidate model worked consistently well, and analysis of the relative importance of different variables revealed that a mix of land use and geographical variables were important independent predictors. Furthermore, several recreation services, outdoor recreation and summer cottages, were poorly predicted by all models. Also, the analysis of model residuals showed that significant spatial correlation remains, demonstrating that important variables are missing. Some part of the unexplained variance is due to things like the geology of the Norrstrom, which was not captured in the drivers used for this study. For example, the proximity of certain municipalities to the Baltic and major rivers results in rapid nutrient transportation from inland to the Baltic Sea, not allowing for high nutrient retention. In this case, biophysical and topographical characteristics of the landscape are key elements Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ for understanding the interactions between services. However, we are less certain about what driver variables would be needed to better predict certain services, like recreation. Most cultural services are the result of the direct experience of nature by humans and can therefore depend on more intangible factors such as individual human preferences (Martin-Lopez et al. 2012, Daniel et al. 2012). Therefore, we suggest that there is a need for eco system service research to develop better integrated social-ecological conceptual models that explain the distribution of ecosystem services in landscapes, and that such theory creation could be informed by theory from geography and other social sciences, especially those that analyze recreation and other cultural uses of landscapes and ecosystems. The distribution of individual services and bundles of ecosystem services across the Norrstrom basin was strongly linked to land use and urbanization (Fig. 3a-c). This region is one of the fastest growing urban areas in Europe (Sporrong 2008), and that explains the importance of Stockholm as a determinant of ES distribution across the study region. One of the striking features of ES in this region is the high level of many services in and near the densely populated city of Stockholm. This might be related with something intrinsic to Swedish culture, where the connection and access to nature assumes a central role for most Swedes (Gidlof-Gunnarsson and Ohrstrom 2007) or some other factors, but our knowledge of what drives the production of cultural services in urban areas is so far poor. A better understanding on what drives these patterns would be useful for future planning in the region, especially identifying how to maintain and enhance these services during its rapid population and urban growth. Our results suggest that land use, and consequently land conversion, and the regional connections to Stockholm are key forces shaping the distribution of ecosystem services. However, there are substantial unknowns about what is determining the distribution of recreation activities and that deserves further investigation. We hope that these results can be developed in further studies to provide some tractable tools to aid in understanding, discussing, and planning for ecosystem services in this rapidly urbanizing region. We expect that one of the most useful ways of achieving this task will be to compare Stockholm region to other regions in Sweden and the world. The approach presented here bridges modeling of ecosystem services that is based on theoretical assumptions, but may not be empirically based, with mapping that is not theoretically based and is not included in modeling (see Table A6.1). Comparison among multiple regions is essential to develop a richer, generalizable understanding of the distribution of multiple ecosystem services across landscapes. We argue that because multiple social-ecological factors appear to codetermine the distribution of ecosystem services, there is a need to develop an explicitly social-ecological theory of what predicts the pattern of ecosystem services. This theory should be developed and tested in multiple locations to identify what aspects of ecosystem services are unique and idiosyncratic and which relationships are general. It would also be useful to consider if analyses using indicators that are not spatially linked would reflect similar results. It would be possible to test provisioning services, for example, economically. Furthermore, to understand how social and ecological changes shape the availability of ecosystem services it is necessary to test how the supply of ecosystem services changes over time. We argue that the empirical analysis of patterns is an important complement to the development of process-based models and expert knowledge-driven models of ecosystem services. We propose that identifying better theories and variables to predict the supply of regulating and cultural services, and applying and extending the statistical approaches presented in this paper to other regions in which multiple ecosystem services have been assessed, would be two useful and productive next steps to understand what types of factors predict the distribution of multiple ecosystem services. Responses to this article can be read online at: http://www.ecologyandsociety.org/issues/responses. php/8077 Acknowledgments: This research was conducted under the Social-ecological dynamics of ecosystem services in the Norrstrom basin (SEEN) project, financed by the Swedish Research Council Formas (project number 2012-1058) as well as the Strategic Research Program EkoKlim at Stockholm University. This research contributes to the Program on Ecosystem Change and Society (www.pecs-science.org). LITERATURE CITED Akaike, H. 1974. A new look at the statistical model identification. Automatic Control, IEEE Transactions 19(6):716-723. http://dx. doi.org/10.1109/TAC.1974.1100705 Alonso, W. 1964. Location and land use. Toward a general theory of land rent. Harvard University Press, Cambridge, Massachusetts, USA. http://dx.doi.org/10.4159/harvard.9780674730854 Bennett, E. M., G. D. Peterson, and L. J. Gordon. 2009. Understanding relationships among multiple ecosystem services. Ecology Letters 12(12):1394-1404. http://dx.doi.org/10.llll/ J.1461-0248.2009.01387.X Bivand, R. S., E. Pebesma, and V. Gömez-Rubio. 2013. Applied spatial data analysis with R. Second edition. Springer, New York, New York, USA. http://dx.doi.org/10.1007/978-l-4614-7618-4 Briassoulis H. 2000. Analysis of land use change: theoretical and modeling approaches. In S. Loveridge and R. W. Jackson, editors. The web book of regional science. West Virginia University, Morgantown, Virginia, USA. [online] URL: http://www.rri.wvu. edu/WebBook/Briassoulis/contents.htm Burkhard, B., F. Kroll, F. Müller, and W. Windhorst. 2009. Landscapes' capacities to provide ecosystem services—a concept for land-cover based assessments. Landscape Online 15:1-22. Burnham, K. P., and D. R. Anderson. 2004. Multimodel inference understanding AIC and BIC in model selection. Sociological Methods & Research 33(2):261-304. http://dx.doi. org/10.1177/0049124104268644 Carpenter, S. R., H. A. Mooney, J. Agard, D. Capistrano, R. S. DeFries, S. Diaz, T. Dietz, A. K. Duraiappah, A. Oten-Yeboah, H. M. Pereira, C. Perrings, W. V. Reid, J. Sarukhan, R. J. Scholes, and A. Whyte. 2009. Science for managing ecosystem services: Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ beyond the Millennium Ecosystem Assessment. Proceedings of the National Academy of Sciences 106(5):1305-1312. http://dx. doi.org/10.1073/pnas.0808772106 Chan, K. M. A., M. R. Shaw, D. R. Cameron, E. C. Underwood, and G. C. Daily. 2006. Conservation planning for ecosystem services. PLoS Biology 4(ll):e379. http://dx.doi.org/10.1371/ journal.pbio.0040379 Chertow, M. R. 2000. The IPAT equation and its variants. Journal of Industrial Ecology 4(4):13-29. http://dx.doi.org/10.1162/1088-1980052541927 Chevan, A., and M. Sutherland. 1991. Hierarchical partitioning. American Statistician 45:90-96. Costanza, R., R. d'Arge, R. de Groot, S. Farber, M. Grasso, B. Hannon, K. Limburg, S. Naeem, R. V. O'Neill, J. Paruelo, R. G. Raskin, P. Sutton, and M. van den Belt. 1997. The value of the world's ecosystem services and natural capital. Nature 387 (15):253-260. http://dx.doi.org/10.1038/387253a0 Costanza, R., and C. Folke. 1997. Valuing ecosystem services with efficiency, fairness and sustainability as goals. Nature's Services: Societal Dependence on Natural Ecosystems 49-70. Cutler, D. R., T. C. Edwards Jr, K. H. Beard, A. Cutler, K. T. Hess, J. Gibson, and J. J. Lawler. 2007. Random forests for classification in ecology. Ecology 88(11):2783-2792. http://dx.doi. org/10.1890/07-0539.1 Daily, G. C, S. Polasky, J. Goldstein, P. M. Kareiva, H. A. Mooney, L. Pejchar, T. H. Ricketts, J. Salzman, and R. Shallenberger. 2009. Ecosystem services in decision making: time to deliver. Frontiers in Ecology and the Environment 7:21 -28. http:// dx.doi.org/10.1890/080025 Daniel, T. C, A. Muhar, A. Arnberger, O. Aznar, J. W. Boyd, K. M. A. Chan, R. Costanza, T. Elmqvist, C. G. Flint, P. H. Gobster, A. Gret-Regamey, R. Lave, S. Muhar, M. Penker, R. G. Ribe, T. Schauppenlehner, T. Sikor, I. Soloviy, M. Spierenburg, K. Taczanowska, J. Tam, and A. von derDunk. 2012. Contributions of cultural services to the ecosystem services agenda. Proceedings of the National Academy of Sciences 109(23):8812-8819. http:// dx.doi.org/10.1073/pnas.1114773109 de Groot, R. S., R. Alkemade, L. Braat, L. Hein, and L. Willemen. 2010. Challenges in integrating the concept of ecosystem services and values in landscape planning, management and decision making. Ecological Complexity 7(3):260-272. http://dx.doi. org/10.1016/j .ecocom.2009.10.006 Dietz, T, E. A. Rosa, and R. York. 2007. Driving the human ecological footprint. Frontiers in Ecology and the Environment 5 (1):13-18. http://dx.doi.org/10.1890/1540-9295(;2007)5[13:dthef| 2.0.co:2 Dormann, C. F, J. M. McPherson, M. B. Araujo, R. Bivand, J. Bolliger, G. Carl, R. G. Davies, A. Hirzel, W. Jetz, W. D. Kissling, I. Kiihn, R. Ohlemuller, P. R. Peres-Neto, B. Reineking, B. Schroder, F. M. Schurr, and R. Wilson. 2007. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography 30:609-628. http://dx.doi.org/10.1111/ J.2007.0906-7590.05171.X Dray, S., and A. B. Dufour. 2007. The ade4 package: implementing the duality diagram for ecologists. Journal of Statistical Software 22(4):l-20. http://dx.doi.org/10.18637/jss. v022.i04 Gidlof-Gunnarsson, A., and E. Ohrstrom. 2007. Noise and well-being in urban residential environments: the potential role of perceived availability to nearby green areas. Landscape and Urban Planning 83:115-126. http://dx.doi.Org/10.1016/j.landurbplan.2007.03.003 Griffith, D. A. 1999. Statistical and mathematical sources of regional science theory: map pattern analysis as an example. Papers in Regional Science 78(l):21-45. Grossman, G. M., and A. B Krueger. 1995. Economic growth and the environment. Quarterly Journal of Economics 110(2):353-377. http://dx.doi.org/10.2307/2118443 Hanspach, J., T. Hartel, A. I. Milcu, F. Mikulcak, I. Dorresteijn, J. Loos, H. von Wehrden, T. Kuemmerle, D. Abson, A. Kovacs-Hostyanszki, A. Bald, and J. Fischer. 2014. A holistic approach to studying social-ecological systems and its application to southern Transylvania. Ecology and Society 19(4):32. http://dx. doi.org/10.5751/es-06915-190432 Harrison, P., and F. Pearce. 2000. AAAS atlas of population & environment. University of California Press, Berkeley, California, USA. Hersperger, A. M., M. Gennaio, P. H. Verburg, and M. Biirgi. 2010. Linking land change with driving forces and actors: four conceptual models. Ecology and Society 15(4):1. [online] URL: http://www.ecologyandsociety. org/vol 15/iss4/art 1 / Jorgenson, A. K. 2003. Consumption and environmental degradation: a cross-national analysis of the ecological footprint. Social Problems 50:374-394. http://dx.doi.Org/10.1525/sp.2003.50.3.374 Kuznets, S. 1955. Economic growth and income inequality. American Economic Review 45:1-28. Lambin, E. F, H. Geist, and R. R. Rindfuss. 2006. Introduction: local processes with global impacts. Pages 1-8 in E. F. Lambin and H. Geist, editors. Land-use and land-cover change: local processes and global impacts. Springer, Berlin, Germany, http:// dx.doi.org/10.1007/3-540-32202-7 1 Liaw, A., and M. Wiener. 2002. Classification and regression by randomForest. R News 2(3):18-22. MacNally, R. 2002. Multiple regression and inference in ecology and conservation biology: further comments on identifying important predictor variables. Biodiversity and Conservation 11:1397-1401. http://dx.doi.Org/10.1023/A:1016250716679 MacNally, R., and C. J. Walsh. 2004. Hierarchical partitioning public-domain software. Biodiversity and Conservation 13:659-660. http://dx.doi.Org/10.1023/B:BIOC.0000009515.11717.0b Martin-Lopez, B., I. Iniesta-Arandia, M. Garcia-Llorente, I. Palomo, I. Casado-Arzuaga, D. G. Del Amo, E. Gomez-Baggethun, E. Oteros-Rozas, I. Palacios-Agundez, B. Willaarts, J. A. Gonzalez, F. Santos-Martin, M. Onaindia, C. Lopez-Santiago, C. Montes. 2012. Uncovering ecosystem service bundles through social preferences. Plos ONE 7:e38970. http://dx.doi. org/10.1371/journal.pone.0038970 Ecology and Society 21(1): 14 http://www.ecologyandsociety.org/vol21/issl/artl4/ Mol, A. P. 2003. Globalization and environmental reform: the ecological modernization of the global economy. MIT Press, Cambridge, Massachusetts, USA. Nelson, E., G. Mendoza, J. Regetz, S. Polasky, H. Tallis, D. R. Cameron, K. M. A. Chan, G. C. Daily, J. Goldstein, P. M. Kareiva, E. Lonsdorf, R. Naidoo, T. H. Ricketts, and M. R. Shaw. 2009. Modeling multiple ecosystem services, biodiversity conservation, commodity production, and tradeoffs at landscape scales. Frontiers in Ecology and the Environment 7:4-11. http://dx.doi. org/10.1890/080023 Nelson, E., S. Polasky, D. J. Lewis, A. J. Plantinga, E. Lonsdorf, D. White, D. Bael, and J. J. Lawler. 2008. Efficiency of incentives to jointly increase carbon sequestration and species conservation on a landscape. Proceedings of the National Academy of Sciences 105:9471-9476. http://dx.doi.org/10.1073/pnas.0706178105 Prasad, A. M., L. R. Iverson, and A. Liaw. 2006. Newer classification and regression tree techniques: bagging and random forests for ecological prediction. Ecosystems 9(2): 181-199. http:// dx.doi.org/10.1007/sl0021-005-0054-l Putz, F. E., and C. Romero. 2001. Biologists and timber certification. Conservation Biology 15:313-314. http://dx.doi. org/10.1046/j.l523-1739.2001.015002313.x Queiroz, C, M. Meacham, K. Richter, K., A. V. Norström, E. Andersson, J. Norberg, and G. Peterson. 2015. Mapping bundles of ecosystem services reveals distinct types of multifunctionality within a Swedish landscape. Ambio 44(1):89-101. http://dx.doi. org/10.1007/s 13280-014-0601 -0 Raudsepp-Hearne, C, G. D. Peterson, and E. M. Bennett. 2010. Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proceedings of the National Academy of Sciences 107 (ll):5242-5247. http://dx.doi.org/10.1073/pnas.0907284107 R Core Team. 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, [online] URL: https://www.r-project.org/ Reyers, B., R. Biggs, G. S. Cumming, T. Elmqvist, A. P. Hejnowicz, and S. Polasky. 2013. Getting the measure of ecosystem services: a social-ecological approach. Frontiers in Ecology and the Environment ll(5):268-273. http://dx.doi.org/10.1890/120144 Schröter, M., R. P. Remme, and L. Hein. 2012. How and where to map supply and demand of ecosystem services for policy-relevant outcomes? Ecological Indicators 23:220-221. http://dx. doi.org/10.1016/j.ecolind.2012.03.025 Sporrong, U. 2008. Features of Nordic physical landscapes: regional characteristics. Pages 568-585 in K. R. Olwig and M. Jones, editors. Nordic landscapes: region and belonging on the northern edge of Europe. University of Minnesota Press, Minneapolis, Minnesota, USA. Stern, P. C, O. R. Young, andD. Druckman, editors. 1992. Global environmental change: understanding the human dimensions. National Academies Press, Washington, D. C, USA. Tilman, D., K. G. Cassman, P. A. Matson, R. Naylor, and S. Polasky. 2002. Agricultural sustainability and intensive production practices. Nature 418:671-677. http://dx.doi. org/10.1038/nature01014 Turner, K. G, M. V. Odgaard, P. K. Bocher, T. Dalgaard, and J.-C. Svenning. 2014. Bundling ecosystem services in Denmark: trade-offs and synergies in a cultural landscape. Landscape and Urban Planning 125:89-104. http://dx.doi.Org/10.1016/j. landurbplan.2014.02.007 von Thünen, J. R. 1826. The isolated state. English edition 1966. Pergainon, Oxford, UK. Wickham H. 2009. ggplot2: elegant graphics for data analysis. Springer, New York, New York, USA. York, R., E. A. Rosa, and T. Dietz. 2003. Footprints on the earth: the environmental consequences of modernity. American Sociological Review 68:279-300. http://dx.doi.org/10.2307/1519769 Zhang, W, T. H. Ricketts, C. Kremen, K. Carney, and S. M. Swinton. 2007. Ecosystem services and dis-services to agriculture. Ecological Economics 64:253-260. http://dx.doi.Org/10.1016/j. ecolecon.2007.02.024 Appendix 1. Detailed methods for data collection of ecosystem services These ecosystem services were assessed in conjunction with the paper, Queiroz et al, 2015. Wheat production The agricultural area set aside for wheat production in each municipality is the indicator of this service. Data was retrieved from the Swedish Board of Agriculture National Statistical Database (Jordbruksverketa). This data is online and publically available. We compiled data for all 62 municipalities for the most recently available year (2012). The total area of wheat production in each municipality was divided by total municipality land area. Cattle The number of cattle on agricultural holdings per municipality in June 2010 is the indicator. This data was available at the Swedish Board of Agriculture National Statistical Database (Jordbruksverketb). The database holds up-to-date data on all cattle existing in Sweden, as cattle owners are obliged to register all animals at the Swedish Board of Agriculture Central Animal Database (CDB). Agricultural holdings refer to activities such as agriculture, stock-farming or horticulture that are undertaken under a single management. Each agricultural holding fulfilled one of the following criteria: a) a minimum of 2 hectares of cropland b) a minimum of 5 hectares of agricultural land c) engaged in commercial horticulture of at least 2,500 km2 outdoor area d) engaged in commercial horticulture with a greenhouse area of at least 200 km2 or: owned a herd that any time from 1 January 2010 to 10 June 2010 included e) at least 10 cattle or f) at least 10 sows, or g) at least 50 hogs, or h) at least 20 sheep or i) at least 1 000 poultry. The total number of cattle in each municipality was weighted by municipality land area. Sheep and pig production The data represents the number of pigs and sheep on agricultural holdings per municipality in June 2010. This data was available at the Swedish Board of Agriculture National Statistical Database (Jordbruksverketb). For sheep and pigs the statistics were based on a survey made to all agricultural holdings in Sweden in June 2010. Agricultural holdings refer to activities as agriculture, stock-farming or horticulture that are undertaken under a single management. Each agricultural holding fulfilled one of the following criteria: a) a minimum of 2 hectares of cropland b) a minimum of 5 hectares of agricultural land c) engaged in commercial horticulture of at least 2,500 km2 outdoor area d) engaged in commercial horticulture with a greenhouse area of at least 200 km2 or: owned a herd that any time from 1 January 2010 to 10 June 2010 included e) at least 10 cattle or f) at least 10 sows, or g) at least 50 hogs, or h) at least 20 sheep or i) at least 1 000 poultry. The total number of cattle in each municipality was weighted by municipality land area. Forest products Products of commercial interest provided by forests such as timber, fiber or paper pulp are considered forest products. The indicator was the area of commercial forest. This indicator represents the potential for the supply of forest products and not actual production values. Data was collected from an online database of the Swedish Forestry Agency (Swedish Forestry Agency), and represents the area in km2 of forest used for commercial purposes in each municipality for the most available recent year (2012). Production forest area in each municipality was weighted by total land area in the municipality. Crop pollination The amount of habitats preferred by pollinators that, within a buffer zone of 200m, intersect cropland area is the indicator of crop pollination. This selects the areas where both the source habitat and potential target crops co-exist. The distance of 200m has been referred in previous studies as a reasonable assumption of the distance a pollinator can fly between its source habitat and the feeding areas. For identifying cropland areas and pollinator preferred habitats we used the Swedish land cover data 2006 (Lantmateriet), which is based on the classification of the Corine Land Cover 2006, a land cover classification for all Europe. Data was retrieved from a governmental online database in raster format. We selected the habitat classes preferred by pollinators in the Marktackedata based on expert personal communication. The habitat classes classified as pollinator habitat and their respective code is listed as it follows: 1.1.2.1.2 Urban areas with more than 200 inhabitants and large green areas 1.1.2.2 Urban areas with less than 200 inhabitants 1.1.2.3 Rural areas with open land 1.2.2 Road and rail networks surroundings. 1.3.1.1 Gravel and sandpits 1.4.1 Urban green areas 1.4.2 Sport and recreational areas 2.2.2 Agricultural land: cultivated berries and fruit 2.3.1 Grasslands 3.2.1 Natural grasslands 3.2.2 Moors and heathland 3.3.1 Beaches, dunes, sand 3.3.4 Burnt areas Raster files were converted to vector format on Quantum GIS geographical information systems software. We intersected the land cover map obtained for all Sweden with the study area, for obtaining an individual land cover map for the Norrstrom basin. Further, we selected all pollinator habitats and performed a buffer analysis (200m). Finally, we intersected the buffer zones with cropland habitat, obtaining the area in km2 with high crop pollination potential for all municipalities. The total crop pollination area for each municipality was weighted by total land municipality area. Standing and running water quality: Water quality data was retrieved from a publically available database at the Swedish Water Information System (VISS). We used data collected during the period 2004-2008 assessing the quality of surface water for both standing (lakes and coastal areas) and running water (rivers and streams). The data used in this study refers to the ecological status of the water, which evaluates if the quality of the water is suitable to the presence of indicator plant and animal species. The classification of the water quality status is made according with three criteria: 1) biological quality of the water 2) physiological and chemical quality of the water 3) hydro-morphological factors. The biological quality of the water assesses the condition of benthic fauna, fish, macroalgae, macrophytes and plankton communities, through the measurement of specific parameters (detailed information on sampling data collection and the parameters used for assessing the state of each group is available on annexes 1 and 4 of the Swedish Marine and Water Authority's regulations on classification and quality standards regarding surface water (HVMFS a). The assessment of physicochemical water quality includes the measurement of synthetic (or artificial) and non-synthetic (such as metals) substances, acidification and eutrophication levels, light and oxygen conditions (HVMFS b). Note that the assessment of the physicochemical water quality is only conducted in those cases where the biological quality status of the water is classified as "good". Hydro-morphological water quality indicators include continuity, hydrological regime and morphology. These are only assessed if the biological, physiological and chemical quality of the water is classified as "good" (HVMFS c). The combination of these three main criteria originates a classification in 5 levels for ecological quality of surface water: 1. High (all three quality criteria are classified as high status) 2. Good (biological and physicochemical criteria are classified as high status but hydro-morphological criteria have a lower classification) 3. Fair (biological criteria are classified as high status while the other two criteria do not fulfill the requirements for high status classification) 4. Bad (none of the three criteria can be classified as "high status") 5. Very bad (none of the three criteria can be classified as "high status") We used this classification to produce a score for water quality of surface water areas in each municipality ranging from from 1 (very bad status) to 5 (high status). Further, we multiplied the classification given to each water surface by its area, obtaining a weighted water quality classification for each water surface. The final value of water quality attributed to each municipality was the sum of the scores calculated for all water surfaces assessed in that municipality, weighted by the total water surface area in the municipality. This procedure was repeated separately for standing water (lakes, coastal and brackish water) and running water (rivers and streams). P and N retention Retention data for both P and N was retrieved from the online database at the Swedish Environmental Emissions Data (SMED). Nitrogen (N) retention refers to the proportion of nitrogen pollution from farmland and private sewers separated from the gross load by soils, lakes and streams before reaching the sea. Phosphorous (P) retention refers to the proportion of phosphorous from all phosphorus loads separated from the gross load by soils, lakes and streams before reaching the sea. P and N retention were expressed by the following expression: Net load (reaching the sea) = (1 - retention) * Gross load. The total P and N retention value obtained for each municipality was weighted by total land municipality area. Moose hunting Moose hunting data was retrieved from a database at the County Administrative Board for the most recent available year (2012/2013) and consisted in a) moose hunting zones that existed in our study area in GIS format (County administrative board a) and b) number of shootings per moose hunting area in Excel format (County administrative board b). All data was processed, sorted by hunting area and converted to csv format in order to match the identification numbers for the hunting areas. Further, we made a join between the shooting data and the hunting areas in QGIS software and all the hunting areas of the study region were merged together. We rasterized the shooting data using raster size 3000 pixels and assessed that the size was detailed enough to minimize the errors. Further, mean number of shootings per municipality were calculated by using zonal statistics. This value was further divided by total land area for each municipality. Summer cottages Data on summer cottages areas was retrieved for the most recently available year (2010) from the Swedish National Statistics Database (SCB). Data is built on a survey done every five years by the administrative agency for statistics in Sweden (SCB). Summer cottages are classified in the database as areas that consist of at least 50 holiday homes with a maximum of 150 meters between them. Holiday homes considered for this survey are buildings that have taxation code 211, 213 or 221 in the Real Estate Taxation Register: 211 - one or two dwelling unit, undeveloped land for holiday homes 213 - one or two dwelling unit, building value less than SEK 50000 221 - one or two dwelling unit, holiday home The total summer cottage area in each municipality was weighted by total land municipality area. Horseback riding We used number of horses in agricultural holdings for the most recent available year (2010) as an indicator of horseback riding. This data was available at the Swedish Board of Agriculture National Statistical Database (Jordbruksverket c). The data is primarily based on a survey conducted in 2010 to agricultural property owners and the numbers correspond to the exact amount of horses hold by agricultural property on the 10th of June 2010. A targeted survey was also conducted at all Swedish riding schools or "similar organizations/companies that offer horseback riding to the public". In addition, data on number of horses and horse holders was collected from animal protection registry (DSK) for 17 municipalities in Stockholm, Gothenburg and Malmo. The response rate was 86% and, theoretically, the populations should not overlap in the results Agricultural holdings refer to activities such as agriculture, stock-farming or horticulture that are undertaken under a single management. Each agricultural holding fulfilled one of the following criteria: a) a minimum of 2 hectares of cropland b) a minimum of 5 hectares of agricultural land c) engaged in commercial horticulture of at least 2,500 km2 outdoor area d) engaged in commercial horticulture with a greenhouse area of at least 200 km2 or: owned a herd that any time from 1 January 2010 to 10 June 2010 included e) at least 10 cattle or f) at least 10 sows, or g) at least 50 hogs, or h) at least 20 sheep or i) at least 1 000 poultry. The total number of horses in agricultural holdings in each municipality was weighted by total land area of the municipality. Outdoor recreation Outdoor recreation data was based on a dataset containing areas of national interest for outdoor recreation. Data was collected in GIS format (County administrative board a) for the counties covering the study area. We used areas of outdoor recreation described under chapter 3, section 6 and chapter 4, section 2 of the Swedish environmental law, Miljobalken (Swedish environmental law). The data was processed in QGIS 2.0.1-Dufour Geographical Information Systems, and the outdoor recreation areas were intersected by municipality borders in order to get a single value for area of outdoor recreation per municipality. The total area of outdoor recreation per municipality was then weighted by total (land and water) municipality area, as outdoor recreation is not only associated to terrestrial areas but also inland water and coastal areas. Proximity to water or coastline generally substantially increases outdoor recreation value of a certain area. Cross country skiing: Cross country skiing as a recreation service was assessed through data from a Swedish online public database for cross-country ski tracks by administrative unit (Skidsar). We retrieved the number of stations with prepared ski tracks for each of the 62 municipalities. The total number of ski stations reported for each municipality was weighted by municipality land area. Biodiversity appreciation Data for biodiversity appreciation was collected from the species database (Artaportalen), that compiles data on species observation across the country. We used the total number of species observations reported to the database for each one of the 62 municipalities for all years until February 2014 for the following groups: vascular plants, mosses, mushrooms, birds, mammals, amphibians, fish and algae. The number of species reported for each group was summed to obtain a single value for each municipality. Further, we divided the number of reported species for each municipality by the total municipality area. Data quality considerations: The Swedish Artportalen is an open database where any individual can report a species observation. This means that many observations are reported by amateurs, i.e. not following a standard or systematic sampling scheme. Therefore, the number of observations reported within a certain geographical area is not necessarily representative of the actual diversity of species of the area, rather reflecting the level of interest that people living or using the area have on biodiversity. References Artaportalen, URL: http://www.artportalen.se County administrative board a, Accessed December 2013, URL: www.gis.lst.se County administrative board b, 'hunting zones,' Accessed December 2013, URL: www.algdata.se HVMFS a, Swedish Marine and Water Authority, annexes 1 and 4, Accessed September 2013, URL: https://www.havochvatten.se/hav/vagledning-lagar/foreskrifter/hvmfs/hvmfs-201319.html HVMFS b, Swedish Marine and Water Authority, annexes 2, 5, and 6, URL: https://www.havochvatten.se/hav/vagledning-lagar/foreskrifter/hvmfs/hvmfs- 201319.html. HVMFS c, Swedish Marine and Water Authority, annex 3, URL: https://www.havochvatten.se/hav/vagledning-lagar/foreskrifter/hvmfs/hvmfs- 201319.html Jordbruksverketa, 'Markanvanding/ Accessed October 2013, URL:http://sta tistik.sjv.se/Database/Jordbruksverket/Markanvandning/J010S MANV_2012_beskr.pdf Jordbruksverketb, 'Husdjur,' Accessed October 2013, URL:http://statistik.sjv.se/PXWeb/Selection.aspx?px_tableid=JO0103G6.px&px_ path=Jordbruksverkets%20statistikdatabas_Husdjur_Antal%20husdjur&px_la nguage=sv&px_db=Jordbruksverkets%20statistikdatabas&rxid=5adf4929-f548-4f27-9bc9 Jordbruksverket c, 'Databasetree,' Accessed October 2013, URL: http://statistik.sjv.se/Database/Jordbruksverket/databasetree.asp Lantmateriet, 'Marktackedata,' Accessed September 2013, URL: https://www.lantmateriet.se/sv/Kartor-ochgeografisk- information/Kartor/Geografiska-teman/GSD-Marktackedata/ Queiroz, C, Meacham, M., Richter, K., Norstrom, A. V., Andersson, E., Norberg, J., & Peterson, G. (2015). Mapping bundles of ecosystem services reveals distinct types of multifunctionality within a Swedish landscape. Ambio 44(1), 89-101. SCB, Swedish National Statistics Database, 'Fritidshusomraden,' Accessed September 2013, URL: http://www.scb.se/sv_/Hitta-statistik/Statistik-efteramne/Miljo/Markanvandning/Fritidshusomraden/12934/2010A01 /Fritids husomraden-2010-lansvis/ Skidspar, Accessed October 2013, URL: http://www.skidspar.se SMED, Swedish Environmental Emission Data, 'Vatten,' Accessed October 2013, URL: http://www.smed.se/vatten/data/plc5 Swedish environmental law, 'Miljöbalken', Accessed September 2013, URL: http://www.riksdagen.se/sv/DokumentLagar/Lagar/Svenskforfattningssamling /_sfs-1998-808/ Swedish Forest Agency, Accessed October 2013, URL:http://www.skogsstyrelsen.se/en/AUTHORITY/Statistics/ VISS, Swedish Water Information System, Accessed September 2013, URL: http://www.viss.lansstyrelsen.se/Exports.aspx?pluginType=StatisticExport&plu ginGuid=01E8EEA2-8F6F-4424-B195-74B949C47238 Appendix 2. Detailed methods for data collection of drivers Education We used the percentage of the adult population in each municipality that has some amount of secondary education as a proxy for education. This data is online and publically available. We compiled data for all 62 municipalities for the most recently available year (2013) (SCBa). Income We used the median net annual income of people 20 years or older for each municipality. This data is online and publicly available. While income is log normally distributed, the average income is fairly equal across the region and more normally distributed. We compiled data for all 62 municipalities for the most recently available year (2012) (SCBb). Wealth Wealth is calculated as the difference between the total amount of assets, both financial and real, and the total amount of debts. This measure indicates a person's net worth and is part of one proxy for development. It is calculated per person and then an average was calculated for each of the 62 municipalities for the most recently available year (2007) (SCBc). Population density We used the log of population density per square kilometer. We log transformed this data due to the huge range of population densities present between urban and rural areas. This data is online and publicly available. We compiled data for all 62 municipalities for the most recently available year (2013) (SCBd). Arable land Arable land was calculated by taking the percentage of total land area of each municipality that is classified as arable land by the statistic bureau of Sweden. The categories of land use are derived from the real estate assessment register. This data is online and publically available. We compiled data for all 62 municipalities for the most recently available year (2010) (SCBe). Grazing land Grazing land was calculated by taking the percentage of total land area of each municipality that is classified as grazing land by the statistic bureau of Sweden. The categories of land use are derived from the real estate assessment register. This data is online and publically available. We compiled data for all 62 municipalities for the most recently available year (2010) (SCBe). Forest area Forest area was calculated by taking the percentage of total land area of each municipality that is classified as forest by the database, Svenska Marktackedata. The data is online and was at the time accessed through Lantmateriet, now available from Naturvardsverket. We compiled data for all 62 municipalities for the most recently available year (2012) (NVV 2015). Built up land Built up land was calculated by taking the percentage of total land area of each municipality that is classified as built up land by the statistic bureau of Sweden. The categories of land use are derived from the real estate assessment register. This data is online and publically available. We compiled data for all 62 municipalities for the most recently available year (2010) (SCBe). Distance from Stockholm We calculated distance from Stockholm using QGIS and centroids of the municipalities. Slope Slope was calculated using QGIS and 50m altitude raster files supplied by the Swedish University of Agricultural Sciences. This data is available to academic institutions (Lantmateriet 2014). References Lantmateriet, 'Hojddata,' accessed June 2014, http://www.lantmateriet.se/sv/Kartor-och-geograiisk-information/Hojddata/GSD-Hojddata-grid-50-/ NVV 2015, Naturvardsverket, 'Svenska Marktackedata,' accessed June 2014, http: //www, milj odataportalen. naturvards verket. se SCBa, Statistics Sweden (SCB), 'Educational attainment of the population,' accessed June 2014, http://scb.se/en_/Finding-statistics/Statistics-by-subject-area/Education-and-research/Education-of-the-population/Educational-attainment-of-the-population/#c_undefined SCBb, Statistics Sweden (SCB), 'Household budget survey (HBS),' accessed June 2014, www. scb. se/he0201 SCBc, Statistics Sweden (SCB), 'Household finances,' accessed June 2014, http://www.statistikdatabasen.scb.se/pxweb/sv/ssd/START HE HE0104/TillgOversiktReg /?rxid=5f366047-908c-4079-8102-cc837407f34f SCBd, Statistics Sweden (SCB), 'Population density,' accessed June 2014, http://www.statistikdatabasen.scb.se/pxweb/en/ssd/START BE BE0101 BEOIOIC/Beft athetkvkmT/?rxid=004b0773-81dd-46fa-8d88-60a69112cfe5 SCBe, Statistics Sweden (SCB), 'Land use,' accessed June 2014, http://www.statistikdatabasen.scb.se/pxweb/en/ssd/START MI MI0803 MI0803A/Mark anvKn/?rxid=e598092a-3d24-4be6-b3cf-8113d5c06809 Appendix 3. Correlation analysis of pairs of drivers Fig. A3.1 20 40 60 80 10 30 50 15 25 35 12.5 13.5 14.5 0 10 30 ..... i i i i i i ...... ..... Distance _ | 0.06 -o.er -0.7T -0.7T -0.7T -o.er -0.33* -0.12 Is | Forest | 0.24" -0.7T -o.sr -0.7$* -o.ee* -0.6T -0.29* -0.30* • '*- ia*.' ■ •; ■" Slope 0.19 0.10 | 1 0.08 | | 0.03 0.01 -0.5*f* -o.er o z *'*■'».. iji-i,. . Built up | 0.8?*| I o.7tr o./tr* -0.10 -0.3f E s • I/"" '' ' Pop. density j row*] | 0.77*1 I o.6tr* 0.07 -0.07 „ E g>- .s^A'?-'-'" | Education | | Q,8T*| | 0.79*1 0.10 -0.09 if'.' ■ ■' '" 1 \ j Income | 0.11 0.01 12.5 145 I_l_t • •'■ •' ■vSfe'*.'.-"1 '1 1 ■ •••• Wealth 0.18 -0.04 . .... I.&V:' I Grazing | 0.5T 3; • '■• I | i' 1 f- ., jjii'-Q ,.;J 1 ■• • Arable 0 100 200 300 1.0 2.0 3.0 4.0 2 4 6 8 180 220 260 300 0 1 2 3 4 The results from the correlation analysis of pairs of drivers with the Person correlation test. Blue represents positive correlations and red represents negative correlations. All significant relationships have *, **, or *** representing their level of significance. Appendix 4. Random forest ranking of all drivers Fig. A4.1 all accur= 0.66 arablejan Forest built_up distancekm log_pop_den log_wealth grazing income education slope T 3 Ranking of all drivers based on the random forest analysis. Appendix 5. Hierarchical partitioning analysis Fig. A5.1 Color Key Log Population Density Log Wealth Income Education Arable Land Grazing Land Built up Land Forest Land Distance from Stockholm The hierarchical partitioning shows the importance of each driver (y-axis) for predicting each ecosystem service (x-axis). Grey boxes indicate no significance. Colored boxes indicate a significant relationship. Darker colors indicate a high correlation. Appendix 6. Analyses performed and their purpose Table A6.1 Analysis Figure Purpose Correlation analysis Fig. A3.1 To explore correlations between pairs of drivers Multiple linear regression Fig. 4 To compare how well the models predict individual ecosystem services Hierarchical partitioning Fig. A5.1 To identify the relative importance of each driver variable Principal component analysis of instrumental variables (PCAIV) Fig. 3 To determine how the assessed drivers and the assessed ecosystem services are related to one another Random forest Fig. 5 and Fig. A3.1 To compare how each model predicts individual ecosystem services and bundles of ecosystem services