Social-ecological drivers of multiple ecosystem services : what variables explain patterns of ecosystem services across the Norrström drainage basin ?

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 socialecological models of human impact on ecosystems explain patterns of 16 ecosystem services (ES) across the 62 municipalities of the Norrström 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.


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 andRomero 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 humandominated 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 Norrström 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 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 Thünen 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.

Study area and ecosystem services
Our study region is the Norrström drainage basin, which is situated in southcentral Sweden.It covers 22,650 km² and includes two of Sweden's largest lakes, Lake Hjälmaren and Lake Mälaren 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 Mälaren 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 Mälaren, 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 Norrström 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 Norrström 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 humanenvironmental 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.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 andKrueger 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 Thünen 1826, Alonso 1964, Griffith 1999).http://www.ecologyandsociety.org/vol21/iss1/art14/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.

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 R² 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 andAnderson 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-environmentalvariables table and a sites-by-species table .Here we adapted it to analyze the relationship between sites-by-drivers and sites-byecosystem 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 andSutherland 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 R², 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 R².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).

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 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 Mälaren and Hjälmaren 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 Mälaren and Hjälmaren, 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 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 Mälaren and Hjälmaren (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.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 bundles across the basin (R² = 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 croplivestock and Mosaic crop-horses, which included most municipalities around lakes Mälaren and Hjälmaren.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 (R² = 0.58) followed by Location Theory (R² = 0.56).These both performed better than the reference model (R² = 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 Norrström 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 andFolke 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 Norrström, 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 http://www.ecologyandsociety.org/vol21/iss1/art14/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-López et al. 2012, Daniel et al. 2012).Therefore, we suggest that there is a need for ecosystem 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 Norrström 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 (Gidlöf-Gunnarsson and Öhrström 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 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 (Jordbruksverket a).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 (Jordbruksverket b).
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 (Jordbruksverket b).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 (Lantmäteriet) 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 Norrström 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 Malmö.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, Miljöbalken (Swedish environmental law).
The data was processed in QGIS 2.0.1--DufourGeographical 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 (Skidsår).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.

Fig. 1 .
Fig. 1.Bundles of ecosystem services in the Norrström 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 socialecological context and all municipalities belonging to a particular bundle are shaded.

Fig. 2 .
Fig. 2. Distribution of drivers across the Norrström 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.

Fig. 3 .
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 Norrström drainage basin, (c) Distribution of the ecosystem service bundles across these municipalities.

Fig. 4 .Fig. 5 .
Fig. 4. Best models for predicting individual ecosystem services.Each bar graph represents each model's R² value for predicting each ecosystem service.

Table 1 .
Models used to estimate ecosystem services and the specific drivers and proxies used in each model.

Table 2 .
Ecosystem services by category and the units in which they were measured.
, 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 Märktäckedata 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.2Urban areas with more than 200 inhabitants and large green areas