Regional Patterns of Ecosystem Services in Cultural Landscapes

European agricultural landscapes have been shaped by humans to produce marketable private goods such as food, feed, fiber and timber. Land-use intensification to increase provisioning services in such productive landscapes alters the capacity of ecosystems to supply other services (often public goods and services) that are also vital for human wellbeing. However, the interactions, synergies and trade-offs among ecosystem services are poorly understood. We assessed the spatial distribution of the services carbon storage, sediment regulation, water yield, crop production, timber supply, and outdoor recreation in the counties Wetterau and Vogelsberg (Hesse, Germany). These counties represent a gradient from intensive arable land use to more extensive mixed land use systems with domination of grassland and forests. Spatially explicit models were used to map the location and quantity of service supply. We addressed the following questions: (1) Where are areas of high and low supply of individual and multiple ecosystem services? (2) Where do the strongest trade-offs and synergies between different services occur? Our results show a pronounced spatial aggregation of different ecosystem services, with locations where at least four services are being supplied at high levels occupying only 5% of the landscape. Indicators for water provision, timber supply, carbon storage, erosion control, and outdoor recreation are positively related to each other, but this relationship is influenced by the trade-offs associated with the ecosystem service food production. Optimization of ecosystem services at the landscape scale has to take these patterns into account.


Introduction
Cultural landscapes have been used by humans for centuries to produce food, fiber, timber and other natural goods, mostly representing private goods. These activities have shaped cultural landscapes by creating specific regional patterns of landscape elements and land use types which reflect the historical and cultural background [1,2]. Alterations of landscapes to maximize production, however, can impact the capacity of ecosystems to provide other-predominantly public-goods and services that are vital for human wellbeing. These include regulating and supporting services, services are distributed across the landscape to determine where interactions such as trade-offs and synergies among ecosystem services occur. We identified hotspots (areas where at least four services were supplied in the upper 25th percentile) and coldspots (areas where at least four services were supplied in the lower 25th percentile) of ecosystem service provision, respectively. We discuss the potential and the limitations of our approach and conclude with recommendations for its application in land use planning and management.

Study Region
The case study area encompasses the counties Wetterau and Vogelsberg located in Hesse, Germany ( Figure 1) and covers an area of 2557 km 2 subdivided into 44 municipalities with areas ranging from 12 to 141 km 2 . The counties represent a typical mosaic of cultural landscapes in central Europe with different climate, topography, soils, and land use.
The county Wetterau (Figure 1, annual precipitation 635 mm, annual average temperature 8.9˝C at Bad Nauheim 148 m a.s.l., period  lies in the rain shadow of the low mountain range Taunus and is distinctly drier than the Vogelsberg [40]. Due to the very fertile loess-rich soils of the Wetterau, the area is predominantly used for intensive agriculture, especially crop production, which comprises more than 40% of the total area, the share of forest and grassland is comparatively small (28% and 14% respectively). This region is one of the oldest cultural landscapes in Germany with evidence for agricultural land use since the Neolithic. Due to its location close to metropolitan area of Frankfurt the county is an economically prospering region with growing industries and a high population density of 268 inhabitants per km 2 [41].

Study Region
The case study area encompasses the counties Wetterau and Vogelsberg located in Hesse, Germany ( Figure 1) and covers an area of 2557 km 2 subdivided into 44 municipalities with areas ranging from 12 to 141 km 2 . The counties represent a typical mosaic of cultural landscapes in central Europe with different climate, topography, soils, and land use.
The county Wetterau (Figure 1, annual precipitation 635 mm, annual average temperature 8.9 °C at Bad Nauheim 148 m a.s.l., period  lies in the rain shadow of the low mountain range Taunus and is distinctly drier than the Vogelsberg [40]. Due to the very fertile loess-rich soils of the Wetterau, the area is predominantly used for intensive agriculture, especially crop production, which comprises more than 40% of the total area, the share of forest and grassland is comparatively small (28% and 14% respectively). This region is one of the oldest cultural landscapes in Germany with evidence for agricultural land use since the Neolithic. Due to its location close to metropolitan area of Frankfurt the county is an economically prospering region with growing industries and a high population density of 268 inhabitants per km 2 [41].
The Vogelsberg (Figure 1) region is a low mountain range of volcanic origin. The higher areas within the Vogelsberg have high precipitation rates (annual precipitation 1307 mm, annual mean temperature 5.6 °C at Hoherodskopf 743 m a.s.l., period 1961-1990) and a relatively long period of snow cover [40]. Due to its geographical characteristics, the agricultural areas in the Vogelsberg are dominated by high share of grassland used for livestock production. The Vogelsberg region has a high share of forest (around 40%), and it has the lowest population density of Hesse with 73 inhabitants per km 2 [41]. Furthermore, the Vogelsberg area is affected by demographic change and high emigration rates to prospering regions such as Frankfurt, leading to local population decline and abandonment of agricultural land. One important source of income is tourism with almost 500,000 overnight stays in the year 2013 [41]. In addition to its role in tourism, this county is one of the largest sources for drinking water in Germany and the main provider of drinking water for Frankfurt and the Rhein-Main area.  [42,43]. Figure 1. Map of the spatial position and land-use of the counties Vogelsberg and Wetterau in Germany including municipality borders [42,43].
The Vogelsberg (Figure 1) region is a low mountain range of volcanic origin. The higher areas within the Vogelsberg have high precipitation rates (annual precipitation 1307 mm, annual mean temperature 5.6˝C at Hoherodskopf 743 m a.s.l., period 1961-1990) and a relatively long period of snow cover [40]. Due to its geographical characteristics, the agricultural areas in the Vogelsberg are dominated by high share of grassland used for livestock production. The Vogelsberg region has a high share of forest (around 40%), and it has the lowest population density of Hesse with 73 inhabitants per km 2 [41]. Furthermore, the Vogelsberg area is affected by demographic change and high emigration rates to prospering regions such as Frankfurt, leading to local population decline and abandonment of agricultural land. One important source of income is tourism with almost 500,000 overnight stays in the year 2013 [41]. In addition to its role in tourism, this county is one of the largest sources for drinking water in Germany and the main provider of drinking water for Frankfurt and the Rhein-Main area.

Assessment of Ecosystem Service Supply
Ecosystem services were classified according to the conceptual framework developed by the initiative 'Common International Classification of Ecosystem Services (CICES)' [44]. The selection of ecosystem services for this study was based on a priority list of ecosystem services in Germany [45] and on the availability of spatial data. In this paper, only ecosystem service supply was considered irrespective of demand [35]. Table 1 provides a summary of the indicators assessed. We used publicly available datasets that can commonly be obtained for this type of landscape to derive generally applicable models. The land cover vector dataset ALKIS served as a basis for the indicators that we used in this study [42]. For the agricultural land cover types a finer classification was possible, overlaying the data with the Integrated Administration and Control System IACS data (InVeKoS, 43). Land cover types were split into 20 classes ( Figure 1, see Appendix A) and converted to a raster dataset of 10 m spatial resolution. All data were imported into ArcGIS 10.2.2 (ESRI) for visualization and data processing. Statistical analysis were performed using the R statistical software [46] (see Section 2.3). Ecosystem services were assessed at a 10 m spatial resolution using data for 2011 or as close as possible to this date.
The InVEST (Integrated Valuation of Ecosystem Services and Trade-offs) modeling tool [47,48] was used to assess the ecosystem services water provision, carbon storage and erosion control, while analyses in ArcGIS 10.2.2 were applied to evaluate food production, timber supply and outdoor recreation. Accuracy of ecosystem service estimates was assessed by comparing our model results with observed measurements (water provision) or census data for this region (timber supply, carbon storage). The ecosystem services food production, erosion control and recreational potential were assessed using regional adapted methods. Full details on methods to derive indicator values are provided in the following paragraphs. A list of spatial datasets used for the assessment of ecosystem service supply is provided in the supplementary material.

Water Provision
We mapped surface water yield as an indicator for the ecosystem service freshwater supply using InVEST. Surface water yield is defined as precipitation minus storage and evapotranspiration losses and therefore describes the amount of water that is discharged as surface flow [47]. The water yield model determines the net hydrological balance as the difference between precipitation and evapotranspiration [47,49], dependent on the characteristics of the vegetation and without distinguishing between surface, subsurface and base flow. In the InVEST 3.0.0 model version we used, the annual water yield pYxjq provided from a cell in the landscape px " 1, 2, . . . , Xq is calculated as follows: where AETxj is the actual annual evapotranspiration on pixel x for landuse and land cover (LULC) j.
Px is the annual precipitation on pixel x. For detailed description of methodology see [47,50]. Annual average precipitation and reference evapotranspiration in the study region were obtained from the German Weather Service [51] for the period 1995-2004. From the soil map (vector dataset) of Hesse [52] root restricting layer depth and plant available water content were derived. The latter was obtained by dividing the fraction of volumetric field capacity by soil depth [47]. Both parameters were converted to a raster dataset of 10 m spatial resolution. For each LULC category specific evapotranspiration coefficients Kc and the root depth were determined based on coefficients reported in the literature (see Appendix A). The empirical constant Z that captures the local precipitation pattern and hydrogeological characteristics was determined according to the following equation [53,54]: where P and AWC are average values of precipitation and available water capacity in the study area. The empirical catchment parameter ω was taken from literature [53]. Based on estimates of our regional catchments [55] we calculated a seasonality factor Z of 4. Validation of the water yield application was performed by comparing InVEST water yield results with aggregated observed runoff values [56]. Observed surface flow for twelve sub-watersheds distributed across and within the study area for the period from 1995 to 2004 was provided by the Hessian Agency for Nature Conservation, Environment and Geology (HLUG). To validate the model, we compared the results through the calculation of the Nash-Sutcliff model efficiency index [57]. The obtained index result was 0.62 (R 2 = 0.70), which indicates that InVEST represented the water flow in the study area sufficiently well.

Timber Supply
A key provisioning service of forests is the supply of timber and fuel wood [58]. Large parts of the study area are covered by forests, and timber production is the most important source of income for forest owners in Hesse [59]. Ownership of forests in Hesse is diverse. In the study area, around 63% of the forest area is managed by the Hessian State Forest Enterprise, while the remaining 37% are privately owned.
The timber stock estimation for the Hessian State Forest Enterprise managed forests were based on national standardized forest inventory data [60]. The inventory provides information on stand area and tree parameters such as species, age, height, timber volume (m 3¨h a´1) and growth rates (m 3¨h a´1¨year´1) on a vector base. The Hessian State Forest Enterprise reclassified forest cover into forest stand types based on the dominating tree species: spruce (Picea), pine (Pinus), Douglas fir (Pseudotsuga), oak (Quercus), beech (Fagus), mix of cherry-maple-ash-elm (Prunus, Acer, Fraxinus, Alnus), mix of birch-willow-poplar-rowan (Betula, Salix, Populus, Sorbus). Information on forest stand-age, coarse wood stocks in cubic meter of wood and average total growth per year were obtained from the forest inventory data. For the privately owned forests that are not managed by the Hessian State Forest Enterprise we estimated the average coarse wood stocks of deciduous (239 m 3¨h a´1), coniferous (312 m 3¨h a´1) and mixed forest (275 m 3¨h a´1) based on the data of the forest inventory. We assigned this information to privately owned forest stands of the LULC data. All timber stock information (including both forest stands managed by the Hessian State Forest Enterprise and privately managed) was converted to a raster dataset of 10 m spatial resolution.
We compared our results to estimates obtained by the Third German National Forest Inventory [61] at federal state level. For this assessment, the average calculated timber stock was calculated over the whole study region. According to the National Forest Inventory [61] the average timber stock in forests in Hesse is 276 m 3¨h a´1; we calculated a similar average of 288.59 m 3¨h a´1 for timber stock in the counties Vogelsberg and Wetterau.

Food Production
To assess the potential of the landscape to produce food we used a frequently applied approach [31,62] of mapping the soil fertility of arable land. Food production depends heavily on natural ecosystem conditions as well as on the management by farmers. For this ecosystem service, the challenge is to separate the relative roles of natural soil conditions and technical measures such as liming, fertilization, plant protection and irrigation. Only the contribution of natural conditions is considered as part of the ecosystem service food production [31,63] and included here. The production of fodder on grassland for animal consumption is not considered in this indicator. To assess the natural potential of the soil to produce crops we mapped potential soil fertility of arable land based on a method recommended for this region of Friedrich and Vorderbrügge [64] using vector based digital soil maps [52] in which water holding capacity, soil moisture and carbonate content were classified into categories and subsequently aggregated to five soil fertility classes. These soil fertility information was subsequently converted to a raster dataset of 10 m spatial resolution. To assess the soil fertility classes of arable fields, we overlaid the soil fertility map with arable land pixels of the LULC dataset.

Carbon Storage
We estimated the amount of carbon stored in the study region by combining the InVEST carbon model with detailed information on wood biomass of forests. The InVEST-toolbox provides a theoretical but comparatively simple production function to assess the capacity of ecosystems to store carbon by summing three major carbon pools: aboveground biomass, belowground biomass and carbon stored in the soil [48]. The carbon pool dead organic matter was not included in this study. The main inputs into the model were a land use map together with a database of the amount of carbon stored for each LULC category for the three carbon pools above-and belowground biomass and soil to 30 cm depth (see Appendix A). These input values were derived from the National Inventory Report for the German Greenhouse Gas Inventory 1990-2011 [65]. Since detailed information on wood biomass of forest stands are available (see Section 2.2.2 Timber Supply) we added this information to the carbon model to derive total carbon stocks across the landscape.
To estimate the carbon stocks in aboveground and belowground biomass of forested areas, the timber volume estimations of national standardized forest inventory [60] was used. For privately managed forests that are not covered by the aforementioned inventory, we calculated average coarse wood stocks, according to the following equation (Equation (3)), following Pistorius et al. [66]: where C is the total carbon in biomass for each forest category; D is the growing stock volume (m 3¨h a´1); rd is the stem bulk density; BEF is the biomass expansion factor; R is the root to shoot ratio and CF is the carbon fraction of dry matter. The factors used to convert raw-wood volume D to dry matter, the bulk density rd, was taken from Knigge and Schulz [67]. Since the available inventory do not contain information about tree biomass other than coarse wood stocks, we applied the literature based biomass expansion factors BEF to account for the biomass of small wood and needles, respectively [65,68]. The root to shoot ratio R for trees was estimated according to the German inventories under UNFCCC [65]. To calculate the carbon fraction CF of biomass, density was set to 0.5 which is used in accordance with the IPCC Guidelines [69]. The carbon stored in forest biomass was added to the above estimates on non-forest biomass and mineral soil carbon stock to provide estimates of the carbon stored across the landscape.

Erosion Control
The ability of the landscape to retain soil is mainly determined by the intensity of rainfall, soil properties, topography, vegetation, and anthropogenic factors such as agricultural practices [70]. Permanent vegetation types act as sinks for sediment including nutrients, pesticides and other agro-chemicals as they interrupt the movement from upslope agricultural areas to surface waters [31,36,71].
Using the Universal Soil Loss Equation (USLE) [70] implemented in the Sediment Delivery Ratio module of InVEST [47] and regional recommendations of input variables, the potential soil loss of each land use cell was calculated as: where A is the potential long-term annual soil loss (t ha´1), R is the erosivity factor of the rain, K is the soil erodibility factor, LS is the slope length and steepness factor, C is the land cover management factor, and P is the supporting practice factor [72]. The value of the retention service is computed by finding the difference between the sediment export from a hypothetical scenario where all land is cleared to bare soil (RKLS) and potential soil loss (A) of the current landscape. To meet regional standardizations we followed guidelines of the Hessian Agency for the Environment and Geology [73] to calculate the factors for the USLE. The rainfall erosivity map was generated for the study region following the method of Schwertmann et al. [74]: where Precip Summer is the long-term (1981 to 2011) average precipitation of the summer months May to October obtained from the German Weather Service [51]. Soil erodibility was estimated by combining characteristics of soil type, organic matter content, and soil texture using the method of DIN 19708 [75]. The slope length and steepness factors LS was calculated using a digital elevation model [55].The land cover specific C factor was assigned using values obtained from the literature (see Appendix A) on a more detailed land use map which discriminates between crop types [43]. To detect areas, which contribute to retention, the Sediment Delivery Ratio model computed the so-called sediment retention index as avoided soil loss by the current land use compared to bare soil.

Outdoor Recreation
Outdoor recreation values are defined as the ability of landscapes to provide opportunities for nature-based recreational activities such as hiking, outdoor sports, camping, bird-watching and fishing [62]. The ecosystem function 'recreational potential' was mapped through a composite indicator based on Paracchini et al. [76]. According to findings from surveys and from literature, recreation potential depends on three different components of people´s behavior and preference for outdoor recreation: the first relates to peoples preference for more natural areas and concerns the degree of naturalness [77]; the second refers to protected areas as they are considered indicators of high natural value [78]; the third component refers to the attractiveness of water bodies [79].
The degree of naturalness is modeled based on the hemeroby index which measures the human influence on landscapes and vegetation and ranges from 1 (natural-without actual human impact) to 7 (artificial) [1]. The naturalness was obtained by attributing to each land cover class its degree of naturalness based on the European hemeroby map of Paracchini and Capitani [77]. Protected areas were mapped using information from the Natura 2000 database of the Birds and the Habitat Directive [80] and other designated areas defined in Germany's Federal Nature Conservation Act [81] such as Landscape Protection Areas, Nature Conservation Areas and Nature Parks [82]. Other designated areas such as National Parks and Biosphere Reserves were not present in the case study area. The attractiveness of water bodies was estimated by calculating the distance to all surface water bodies. We included all stagnant water bodies extracted from the LULC dataset and water courses of the first and second order [42]. Small streams of the third order where not included.
We assumed that each component covers a different aspect of outdoor recreation. After normalizing each variable to a range to 0-1 using a Min-Max normalization method [83] the three components have been aggregated considering each variable equally important therefore given equal weights [76]. We identified areas with high recreational potential, hereinafter "recreational hotspots", by extracting the patches within the upper 25th percentile [18]. This analysis addressed recreation and not tourism, ranging from a short walk or bicycle ride to a Sunday hiking trip. Long distance travelling (>100 km) was not included in the study as the current analysis focused on resident populations.

Spatial Relationships among Ecosystem Services
Individual ecosystem services were mapped to visualize and compare their spatial patterns. We calculated summary statistics and evaluated the degree for spatial clustering of each service using Moran's I of package "raster" [84]. We identified service provisioning hotspots and coldspots of each service and the spatial congruence of those for multiple service supply. Grid cells in the study area that were within the upper 25th percentile were defined as individual hotspots, and those within the lowest 25th percentile were considered as individual coldspots. Hotspots and coldspots of multiple services were defined as those areas where at least four services were supplied in the upper 25th, respective the lower 25th percentile. Spatial patterns of hotspots and coldspots were analyzed by computing the proportion of the study area occupied, the patch density and the mean patch area in R using ClassStats of the package "SDMTools" [85,86].
Analyzing the interactions between pairs of ecosystem services took place on a municipal level. The municipal boundaries were taken from the Administrative Map of Germany [87]. The area of each municipality was calculated based on these boundaries. The administrative boundaries were used to identify trade-offs and synergies spatial variables on the supply of ecosystem services. Municipal level was chosen since these administrative units represent the smallest unit of governance and political decision making (such as landscape and land-use planning). Because the municipalities are not all of the same size (area ranges from 12 to 141 km 2 , mean 58 km 2 ) we aggregated the sets of ecosystem service associated with each municipality (n = 44) by extracting the mean values of individual ecosystem services indicators per municipality. Non-normally distributed indicators were transformed (water provision x 0.2 , food production x 0.3 erosion control x 0.5 , recreational potential x 0.2 ) using the Box-Cox transformation of the "MASS" package [88] before statistical analyses. The transformed data were standardized by ranging the data to a scale from zero to one based on minimum and maximum values to allow for meaningful comparisons of all indicators [89]. Graphical and correlation analysis were performed to examine interactions between pairs of ecosystem services. For correlation analysis, the Pearson parametric correlation test was used to assess correlation coefficients for all pairs of ecosystem services. For graphical analysis we used a bivariate version of the boxplot, the bagplot [90]. We calculated bagplots for all possible pair-wise comparisons using the package "aplpack" [91]. Synergetic relationship between two ecosystem services was illustrated in the bagplot when the plot was oriented from the lower left to upper right and thus covering the negative/negative and positive/positive space. Trade-offs between two ecosystem services are expected if the bagplot is oriented from the upper left to the lower right.

Spatial Patterns of Ecosystem Service Supply
Supply of individual ecosystem services (measured at a resolution of 10 mˆ10 m) varied across the landscape (Figure 2) and showed a pronounced spatial clustering (all ecosystem services Moran's I > 0.60). For descriptive statistics of each service, see Table 2. Within the study area on average 100 t¨ha´1 of carbon are stored in the above-and belowground biomass and the soil. The biggest carbon pools can be found in the deciduous forest stands with carbon stocks above 250 t¨ha´1. The sediment model showed an average potential soil loss of around 200,000 t¨year´1 across the study area, with potential soil loss up to 14 kg¨m´2. However, only 4300 t¨year´1 are exported to streams and rivers, because most of the soil particles are retained by vegetation. Calculating the water provision of the study region revealed an average water yield of 255 mm¨year´1. However, these values vary in space ranging from water yields under 100 mm in the Wetter valley up to over 800 mm in the mountainous region of the Vogelsberg with higher annual rainfall and lower evapotranspiration. The total standing stock of timber had a volume of more than 26 Mio. m 3 across the study area with highest timber stocks in coniferous forest stands of 312 m 3¨h a´1. The total area of arable land cells in the two highest soil fertility classes comprises 138 km 2 which represents around 12% of the counties. Regions of high recreational potential, so-called recreational hotspots, with values within the upper 25th percentile, cover only 1.20% (31 km 2 ) of the study region. All ecosystem services are clustered according to ecological, geographic or social factors that led to land use systems and associated ecosystem services. For example, food production is found in rather flat, fertile regions of the study region, while timber is mainly provided from steeper slopes or higher elevations. There were similarities among the spatial distribution of the supply of different services (e.g., timber supply and carbon storage). However, the geographic distribution of the supply of individual ecosystem service revealed that each pattern was distinct. Timber for example is supplied in forest areas only, while carbon is stored in almost all land use types including agricultural areas. The quantities of supply vary, leading to distinct patterns of the respective ecosystem services. Only 5% of the landscape was occupied by ecosystem service hotspots, where at least four services were provided in the upper 25th percentile (Figure 3a). These areas often coincided with forest areas and other wooded structures in mountainous areas and were few in number and rather large in size (patch density = 7.4 km´2, mean patch area = 60.9 ha). Coldspots, which are areas where at least four services were produced in the lower 25th percentile, occupied around 20% of the landscape (Figure 3b) and corresponded with cropland, urban areas and other sealed areas. Coldspots were more numerous and slightly smaller in size (patch density = 35.3 km´2, mean patch area = 57.3 ha).

Relationships among Ecosystem Services
Analysis of relationships of ecosystem service supply was carried out at the municipal level. Most of the ecosystem services were correlated (Figure 4). Correlations were highly significant for twelve out of 15 pairs of ecosystem services (p-value < 0.001). The relationship between regulating services (carbon storage and erosion control) and the provision of timber was synergistic. The cultural ecosystem service outdoor recreation was found to have a positive relationship with areas of high water supply (r = 0.71). Food production was negatively correlated with all other ecosystem services included in this study, indicating a trade-off. The ecosystem service pairs water provision vs. timber supply (r = 0.25) water provision vs. carbon storage (r = 0.34) and erosion control vs. outdoor recreation (r = 0.37) were not significantly correlated. services (carbon storage and erosion control) and the provision of timber was synergistic. The cultural ecosystem service outdoor recreation was found to have a positive relationship with areas of high water supply (r = 0.71). Food production was negatively correlated with all other ecosystem services included in this study, indicating a trade-off. The ecosystem service pairs water provision vs. timber supply (r = 0.25) water provision vs. carbon storage (r = 0.34) and erosion control vs. outdoor recreation (r = 0.37) were not significantly correlated.

Assessment of Ecosystem Service Supply
We aimed to assess spatial variability in the supply of multiple ecosystem goods and services in the study region. To evaluate ecosystem services in a replicable and standardized manner [92][93][94], we used biophysical proxies based on datasets that are typically available in this type of landscape, and models which can be applied in other locations. The ecosystem service indicators reflect the capacity of ecosystems to provide services based on present knowledge of ecological relationships. However, for many of the ecosystem services we rely on simplified indicators to estimate the natural potential of providing these services [95]. The use of relatively simple ecosystem service models such as InVEST enables the extension to other study locations, the assessment of multiple services and may also generate greater transparency and trust among users [94].
In the models, neither spatial flow dynamics such as down-stream effects of landscape features on ecosystem service performance nor feedbacks between ecosystem services and temporal variability were taken into account. The erosion control model for example does not consider feedbacks of sedimentation processes on water quality and quantity or carbon storage. Moreover, the erosion control model uses long-term average summer precipitation data to derive the erosivity factor. However, extreme precipitation events (either within or outside of summer season) may result in dramatic soil losses. The used model is not able to detect temporal variability of sediment erosion. Further, specific landscape features such as wetlands or hedges can have disproportionate effects on various ecosystem services dependent on their position in the landscape. This information is essential for strategic spatial positioning of key landscape elements to improve multiple ecosystem services [96].
The carbon storage model does not consider the full carbon cycle or fluxes among carbon pools [47]. In addition, it relies on carbon storage estimates derived from land use classes and does not consider other controlling factors for carbon storage such as climate, soil depth or management practices [47,97]. The freshwater model is also based on simplifications of well-known hydrological relationships [27,56] and works on an annual basis. Seasonal variability in the delivery of hydrological ecosystem services is not considered here. Furthermore, we did not include water withdrawal in the freshwater model which could lead to a general inaccurate estimation of water provision [98]. Neglecting ecosystem service flow dynamics and feedbacks may lead to over-or underestimation of ecosystem service supply. Modeling approaches that incorporate service-specific sources, sinks, users, spatial flows and interactions can provide a more complete understanding of ecosystems services [94,99]. Although more complex modeling approaches may more accurately represent natural processes and interactions, the potential gains in accuracy associated with this must be weighed against the increased complexity and reduced generalizability [47].
In applying the models, a conversion of (polygon) vector datasets (land use, soil information, protected areas and forest inventory data, see supplementary material) to raster format was necessary. This transformation, however, may produce considerable distortion in landscape pattern characteristics (e.g., cut cells or loss of small or thin elements) and, hence, in the outcome of the ecosystem service assessment [100]. In order to maintain the integrity of our data and reduce the loss of geometric precision we used a fine spatial resolution of 10 m.
These limitations could lead to inaccurate estimates of ecosystem service supply. For most ecosystem services, however, we were able to assess the accuracy of our estimates by validation (potential water provision) or comparison with other statistical estimations (potential timber supply) or applied the model according to regional knowledge (e.g., potential carbon storage, potential erosion control) to minimize overall uncertainty.

Relationships among Multiple Ecosystem Services
To derive total ecosystem service supply we calculated multiple ecosystem service hotspots and coldspots, respectively, which represents the capacity to provide multiple services. This method does not include the demand for the particular ecosystem service. To improve the method of defining hotspots and coldspots we propose to investigate which ecosystem services are below a critical threshold in relation to demand [35]. Adding up of the supply of ecosystem services to evaluate multifunctionality, does take interaction among ecosystem services indirectly into account (e.g., carbon storage and timber supply) since they are both included independently of each other. Comparison of cumulative ecosystem service supply (represented by hotspots or coldspots) with land use information suggests a strong correlation with forest, woodland and grassland rich regions. The abundance and spatial organization of fields, hedgerows, grassland and forest patches determine to large parts the ecosystem service supply [9]. Our results, however, depend on the set of ecosystem services analyzed. Nonetheless, these results support existing knowledge about links between landscape structure and ecosystem services [17,37,38], proposing that methods applied here can be applied to other cultural landscapes.
Correlations between ecosystem services analyzed in this study were found to be stronger than those between different ecosystem services found in a paper from Jopke et al. [38,101]. This may be due to the much larger European Union wide spatial scale at which this study was conducted (NUTS 3 level across Europe), the number of services assessed and different proxies used. The scale of the analysis will to some extent control the intensity of interactions. We suggest the spatial scale of municipalities to be small enough so that variables that determine the supply of one service will also have similar effects on the supply of other services. The comparison of our outcomes with a study conducted in Quebec, Canada on a municipal level with comparable size of administrative units (averaging 74 km 2 ) revealed similar results [37].
Among all studied ecosystem services pairs, most showed a synergetic relationship, and correlations were significant. Some ecosystem service pairs showed no significant correlation. Erosion control and outdoor recreation showed no significant correlation which may be because the recreational potential of the landscape is determined by the degree of naturalness, designated areas and water areas which are seldom located in arable areas of high erosion control. Water provision was not correlated to timber supply and carbon storage. This may be due to the fact that precipitation patterns, which control water yield, follow a topographic gradient while forest areas with high rates of timber supply and carbon storage are distributed across the study area. All ecosystem services showed trade-offs against crop production. Analyses of interactions among ecosystem services in Quebec, Canada and Europe respectively lead similarly to the conclusion that crop production is the most conflicting ecosystem service [37,38,78]. Our analysis, though, shows only correlational relationships, yet causal relationships cannot be proved [27]. However, there seems to be a general trade-off between crop production and regulating or cultural services. Regions rich in arable production essentially produce crops and are relatively poor in delivering other ecosystem services. Among all other considered services, synergies were found. Regions rich in green infrastructure such as hedges, wetlands and forests provide a wide array of services. Landscape composition and landscape diversity are important determinants of ecosystem services.
The results presented in this manuscript have practical uses for landscape management, spatial planning and policy in the analyzed counties. Spatial explicit assessment of ecosystem service supply will help to identify priority regions for landscape management to promote landscape multifunctionality. Having identified trade-off regions can help decision-making in targeting management actions (e.g., implementation of agri-environmental schemes). To secure the overall supply of ecosystem services the development of green infrastructure should be especially promoted in the intensive agricultural areas of the counties Wetterau and Vogelsberg. Furthermore, the impact on ecosystem provision from a proposed landscape alteration will be more thoroughly evaluated if improvements or declines in analyzed ecosystem services can be demonstrated to stakeholders. Finally, assessment of ecosystem services will be used in reaching decision that consider the full effects of a proposed use of an ecosystem, rather than just those costs or values that enter markets in the form of private goods [96]. Ecosystem service assessment methods like the ones applied should be included in spatial planning in order to identify priority regions for management actions and consider consequences of land use changes for different ecosystem services. Ecosystem service assessment can improve landscape planning and decision-making and control for a sustainable use of ecosystem services.

Conclusions
We provide a methodological framework for mapping and analyzing multiple ecosystem services in a spatially explicit way. Using this framework, we exposed positive, negative and neutral relationships among ecosystem services at the municipal level. Our study demonstrates that ecosystem service supply shows distinct spatial clustering. We found positive relations among water provision, timber supply, carbon storage, erosion control, and outdoor recreation, but trade-offs between food production and all other services. Highlighting these patterns is a prerequisite for optimizing multiple ecosystem services through adaptive land-management decisions. Applying this knowledge in practice requires that social and ecological drivers of change in supply of and demand for ecosystem service are further investigated. A mechanistic understanding of the feedbacks between drivers, ecosystem services and societal responses to ecosystem service change is needed to achieve targeted management schemes that ensure long-term, balanced provision of ecosystem services.
Supplementary Materials: The following are available online at www.mdpi.com/2073-445X/5/2/17/s1. Acknowledgments: This study was supported by the Federal Ministry for Education and Research (BMBF) in the project JAGUAR (Project-ID 01LC1106A). We are grateful to all colleagues working on this project for continuous discussion and support. Finally, we thank three anonymous reviewers, who helped us to improve the quality of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Biophysical table with input parameter to model carbon storage, sediment retention and water yield including References: Carbon stock estimates for above-and belowground biomass and soil to 30 cm depth; land cover management factor C, supporting practice factor P, maximal rooting depth and evapotranspiration coefficients Kc. Table A1. Biophysical table with input parameter to model carbon storage, sediment retention and water yield: Carbon stock estimates for above-and belowground biomass and soil to 30 cm depth; land cover management factor C, supporting practice factor P, maximal rooting depth and evapotranspiration coefficients Kc.