Global change impacts on ecosystem services : a spatially explicit assessment for Europe

Corresponding author: Chiara Polce (chiara.polce@gmail.com) Academic editor: Louise Willemen Received: 23 Jul 2016 | Accepted: 25 Oct 2016 | Published: 1 Nov 2016 Citation: Polce C, Maes J, Brander L, Cescatti A, Baranzelli C, Lavalle C, Zulian G (2016) Global change impacts on ecosystem services: a spatially explicit assessment for Europe. One Ecosystem 1: e9990. https://doi.org/10.3897/oneeco.1.e9990


Introduction
The widely reported impacts of climate change on ecosystems and biodiversity (e.g.Feehan et al. 2009) pose a threat to the supply of ecosystem services (Schröter et al. 2005).Ecosystem services (ES) arise when ecological structures or functions contribute toward meeting a human demand; they are, arguably, underpinned by biodiversity.Global change impacts on biodiversity and ecosystems, therefore, it is likely to also affect the supply of ES and, consequently, the overall human well-being.
Fig. 1 summarises the relationship between ecosystems and human wellbeing (Haines-Young and Potschin 2010, Maes et al. 2013).From left to right, the first two boxes relate solely to the ecosystem dimension and involve the bio-physical processes underpinning ecosystem functioning; for instance, the processes by which atmospheric carbon dioxide is converted into organic carbon compounds such as carbohydrates.The last two boxes relate to the human dimension, and characterise ecosystem functions that satisfy human needs; for instance, the production of aboveground wood biomass that can be sold as timber.
Ecosystem functions are defined as the capacity of natural processes and components to provide goods and services that satisfy human needs, directly or indirectly; accordingly, they form a subset of ecological processes and ecosystem structures (de Groot 1992).Fig. 1 is also known as the 'ecosystem services cascade' since it suggests that ecosystem services flow from ecosystems to society.Especially in Europe, it is a commonly used framework for the assessment of ecosystems and their services.For instance, it is applied by the TEEB (The Economics of Ecosystems and Biodiversity) and by MAES, the EU initiative under Action 5 of the Biodiversity Strategy to map and assess ecosystems and their services.
To date, different classification systems for ES are in use.They invariantly discriminate among provisioning services (the goods we obtain from ecosystems), regulating and maintenance services (the capacity of ecosystems to maintain a liveable environment), and cultural services (the non-material benefits).Well known provisioning services are water, timber, fish and agricultural products which are traded on markets.The supply of provisioning ecosystem services is heavily influenced by human inputs in the form of energy, labour and nutrient subsidies.Regulating services include the removal of pollutants from soil, air and water, or services which support crop production such as pollination and soil erosion control.They are directly dependent on the proper ecological functioning of ecosystems; they are not traded on markets so that their contribution to human well-being is more difficult to assess.Cultural services include, among others, nature-based recreation and tourism; they are, more than other services, essentially defined by human preferences.
A recent assessment of the potential consequences of climate change in Europe on several impact categories (i.e.agriculture, river floods, coastal areas, tourism and human health) has revealed large variations across regions in terms of damages and benefits (Ciscar Martinez et al. 2011).Climate change impacts on ES were only partly addressed by this assessment (they were mainly assessed in relation to agricultural yields) and were identified as one of the next knowledge gap to address.A common approach to assess ES involves the use of land use and land cover (LULC) information in combination with other statistics and data to infer stocks and flows of ES generated by different ecosystems (Maes et al. 2012).LULC data enable us to assess the impacts of changes in land management on ES (Maes et al. 2013); hence they also provide information of policy relevance (e.g.Liquete et al. 2013, Maes et al. 2013, Maes et al. 2012).
Considering both LULC projections and climate change (CC) scenarios, in principle, enables us to capture the main pressures acting on ecosystems and ES, thus enhancing the suitability of this approach to generate policy-relevant information.The establishment of the linkages among CC/LULC/ES, however, requires, amongst others, accurate considerations on issues related to the spatial and temporal scales characterising the involved processes; therefore, it can only be tackled in a stepwise approach -depending on the availability of data, tools and resources.
The focus of this study lies on regulating ecosystem services.This choice follows from differences among the main groups of ecosystem services: the delivery of regulating services, in fact, is tightly coupled to well-functioning ecosystems.This property makes regulating ES suitable examples for a first case study at a Pan-European scale.
This study provides an assessment on the possibility to establish relationships between CC, LULC and changes in ES (potential) supply, via bio-physical parameters that: (i) respond to climatic drivers; (ii) can be tied to ES with semi-quantitative or qualitative relationships, without the need of process-based models (which are beyond the scope of this study).This study is preliminary for a number of reasons: for instance, we only looked at impacts on three ES; additionally, due to the difficulties to assess the role of biodiversity as a service providing unit, we do not look at the impacts of climate and LULC changes on biodiversity, although we acknowledge that biodiversity is affected by changes in these drivers, and that biodiversity underpins ES.With these limitations in mind, our work quantifies, in relative terms, changes in ES provision across terrestrial ecosystems, induced by changes in both LULC and climate and, also, by changes in climate alone.

Overview
Every box in the cascade model of Fig. 1 is dependent on specific expertise and knowledge about the structure of ecosystems (e.g.location, relative share), their functioning (interplay between abiotic and biotic environment, ecological processes), and the services they provide (which requires knowledge about human demand), making ecosystem assessment a complex process.Furthermore, an analysis of benefits and values requires inputs from economists or social scientists.Every step from left to right in the ecosystem services cascade, therefore, involves the integration of different expertise, knowledge, data and models and requires the transfer of ecological values into economic values.To our knowledge, there are no models or expert systems available that can perform these tasks at once, and also in this study we used a pragmatic approach by including different models and linking them to each other.
In addition to selecting regulating services for being tightly coupled to well-functioning ecosystems, the final choice to include certain services was further influenced by the different models that were selected to establish the linkages between CC, LULC and ES:

•
The Community Land Model (CLM) (Oleson et al. 2010) provided the climate based scenarios of changes in ecosystems and ecosystem functioning.CLM examines the physical, chemical, and biological processes by which terrestrial ecosystems affect and are affected by climate.Some of the processes present within CLM are, for instance, vegetation and soil hydrology, energy fluxes, Carbon and Nitrogen cycling, and routing of runoff from rivers to oceans.Thus, CLM was the source of the bio-physical processes and structures as outlined in the cascade model (Fig. 1).

•
The LUISA modelling platform (Baranzelli et al. 2014) was used to spatially identify different ecosystems based on land cover and land use.Thus, LUISA was used to assess the actual ES provision, translating processes and functions into services ( Fig. 1).
The choice for these two models constrained the number of ecosystems and ecosystem services that could be covered in this assessment.As a result, we present the assessment for terrestrial ecosystems based on their capacity to deliver air quality regulation, soil retention (or soil erosion control) and water (quantity) regulation (capacity to regulate water flows which either supply or buffer water).
Fig. 2 is a schematic representation of the main data sources and methods leading to assess ES provision and changes therein.In summary, we selected three regulating ES which we estimated using Bayesian Networks (BN) built with bio-physical components from the CLM.The estimated ES were then spatially linked to the relevant ecosystems, using LUISA projected LULC types for the present baseline ( 2010) and future (2050).After isolating changes due to climate alone, we also highlighted areas with the greatest climateinduced changes in ecosystem services provision, for a selected ecosystem: forest.
Changes in vegetation and their ecological functioning are expected to affect the delivery of ecosystem services.Climax vegetation, in particular forests, are key suppliers of regulating ecosystem services and control the quantities and quality of air, water, soil, and biomass.Hence, changes in the structure of the vegetation (e.g.larger leaf area index), or in the total biomass (e.g. through enhanced growth or fire) are likely to have a detectable effect on the services provided by forests and other ecosystems.Furthermore, forests cover about 42% of the EU-28 mainland: small changes in the supply of ecosystem services, therefore, could result in very high gains or losses in absolute terms.
The next subsections provide a detailed description of the methods we followed.
Global change impacts on ecosystem services: a spatially explicit assessment ...

ES choice and proxies
The ecosystem services (potential) supply was estimated using Bayesian Networks (BN), which enable us to model under uncertainty and to integrate different types of probabilistic information (Cain 2001).These features make BN attractive in environmental management problems, with successful applications to land use decisions, impacts of climate change and assessment of several ES (Barton et al. 2012, Landuyt et al. 2013, Molina et al. 2013, Richards et al. 2013, Celio et al. 2014).The Suppl.material 1 shows a simple BN and it briefly describes its main components; for an introduction to BN for natural resource management, we refer to Cain (2001).
We used expert knowledge and available literature to select outputs from CLM which can be related to regulating ES (i.e. to define the nodes of the BN), as well as to define the relations between the selected CLM components and each specific ES (i.e. to define each ES model).The CLM outputs consisted of projections of terrestrial processes under the IPCC SRES A1B climate scenario, generated by three combinations of global and regional climate models (RCM).The CLM outputs have a spatial resolution of 25 km and a monthly temporal resolution, covering the period 1961 -2099.We selected two 30-year nonoverlapping periods to maximise the temporal overlap between CLM and LUISA frame: the present or 'baseline' (from 1991 to 2020), and the 2050 (from 2021 to 2050).For each period and CLM component, we then derived basic descriptive statistics such as minimum, maximum, average and standard deviation, pooling together the results from the three RCM.This information was used to reclassify the original CLM predicted values into relative scores between 0 and 1, and then to assign them their respective probabilities (see next paragraph).As stated above, the selected CLM components represent different nodes of the BN.Conditional Probability Tables (CPT) were filled to reflect the distribution of states (i.e.their probability, p) resulting from each different climate model.More specifically, for each CLM output of interest, we applied Fisher's intervals (Fisher 1958) between 0 and 1 to classify the full range of values predicted by the three different climate models over the 30-year time frame.We then computed the 30-year average for each RCM, and we assigned it to the corresponding interval.If the predictions from the three models fell all within the same interval, the corresponding state in the CPT was assigned the greatest probability (i.e.p = 1).If there was agreement from two out of the three RCMs, the two corresponding states were chosen, one with p = 0.66, the other one with p = 0.33.
We built BNs for the following regulatory ES: Water flow regulation This allowed us to detect and characterise how processes and functions underpinning the selected ES change over time.In particular, we assessed the projected changes for 2050, in relation to the baseline.
We selected 12 bio-physical variables from CLM, which we used as proxies for ecosystem functions or ecosystem processes, after establishing a positive or negative relationship with each relevant ecosystem service (Table 1).Fig. 3 is a representation of the BN built for the assessment, featuring the linkages between the bio-physical components and the relevant ES.The next section provides a description of the ES models.

ES models
Here we present the ES models used within the BN, with reference to the relevant literature consulted to derive the relationships between the different model components.

Air quality regulation
The BN model to estimate air quality regulation included the contribution of three main elements: vegetation, rainfall and fire.The first two exerted a positive effect on air quality, whilst the last one a negative effect.
• Vegetation: Vegetation leaves absorb gaseous pollutants through their stomata, while particles are removed from the air by deposition onto leaves and branches (Vos et al. 2013).Pollution removal is known to be positively related to tree cover and length of in-leaf season, and to vary with meteorological variables that affect tree transpiration and deposition velocities (Nowak et al. 2006).The contribution of vegetation was rendered by an additive effect of foliage and canopy height (ELAI and HTOP respectively, Table 1).See Table 1 for their definition and see the main text for a description of each relevant ES.
[Eq. 1] • Rainfall: Positive effects on removal capacity are also exerted by rainfall through wet deposition (Cooter et al. 2013), although increased precipitation may also lead to reduced total removal via dry deposition (Nowak et al. 2006).In absence of information to quantify this effect, we accounted for this potential reduction by assigning a weight = 0.8 to the positive contribution of rainfall (RAIN, Table 1).The implication of this choice are addressed in the 'Discussion'.

•
Fire: Wild or prescribed, fire can significantly degrade air quality by impairing visibility and releasing toxic by-products into the atmosphere.The susceptibility to open biomass burning events is likely to increase in Southern Europe, due to warmer and drier conditions predicted by climate change scenarios (Giannakopoulos et al. 2009).Fire-related impacts and potential mitigation options are currently being investigated by several authors (e.g.Garcia-Hurtado et al. 2014, Sarigiannis et al. 2014).We accounted for the negative effect exerted by the fire by subtracting fire probability (FIRE_P, Table 1) to the benefits derived from vegetation and rainfall.
In summary, the BN model for air quality regulation included five parent nodes (Fig. 3), with final scores potentially ranging between -1 and 2.8, where negative values indicate air pollution by fire.These scores were then rescaled to the 0 -1 interval (with the original 0 placed at 0.26, and values < 0.26 indicating air pollution by fire).The relative contribution of each LULC type to regulating air quality is described in the section 'LULC contribution to ES delivery' and in Suppl.material 2. Equation 1 summarises the model:

Soil erosion control
The BN model to estimate soil erosion control included the protective role of vegetation and soil organic content, and the negative effect of surface runoff, which causes erosion by water.We did not consider the potential contribution of specific support practices aimed at controlling erosion, such as policy implementation or voluntary measures, as this would require information we currently do not have.When we scored the contribution of different LULC types for their capacity to control soil erosion, however, we took into consideration that certain man-made elements can potentially provide a positive contribution to controlling erosion (Burkhard et al. 2014).See section 'LULC contribution to ES delivery' for further details on the relative contribution of LULC classes to soil erosion control.
• Vegetation: We chose the amount of Carbon in fine roots (FROOTC, Table 1) to represent the protective role offered by vegetation.This choice was motivated by recent experiments investigating the role of root traits to reduce erosion rates by concentrated flow (Burylo et al. 2012).The results highlighted that the potential to reduce erosion was negatively correlated to root diameter, and positively correlated to the percentage of fine roots.• Soil organic matter: In addition to fine roots, we also included soil organic matter (TOTSOMC, Table 1) as a protective component; greater organic matter into the soil, in fact, provides better structural and water-holding qualities (Lal 1994).

•
Surface runoff: We used surface runoff (QOVER, Table 1) to represent the erosion caused by water flow (Toy et al. 2002).
The estimated soil erosion control resulted from subtracting the effect of surface runoff to the estimated protection obtained from fine roots and soil organic matter.
In summary, the BN model for soil erosion control included four parent nodes (Fig. 3) and rendered scores between -1.0 and +2.0, with negative values indicating erosion and positive values protection from erosion.These scores were then rescaled to the 0 -1 interval (with the original 0 placed at 0.33, and values < 0.33 indicating erosion).Equation 2summarises the model: Where: • Protection = Vegetation (FROOTC)+ Soil organic matter (TOTSOMC) Water flow regulation The BN model to estimate water flow regulation included the contribution of four components of the hydrological cycle, making up the surface and groundwater: rainfall, water loss (from evaporation and transpiration), water in the unconfined aquifer and aquifer recharge rate.We chose these variables as they might be directly affected by changing climate, for instance as a result of changes in temperature or rainfall patterns (Wu et al. 2012).
• Surface water: We chose precipitation amount (RAIN, Table 1) minus water loss, to represent the potential source of surface water (Bagstad et al. 2011).We represented water loss by summing the contribution of ground evaporation (QSOIL, Table 1) and evapotranspiration (canopy evaporation and canopy transpiration, QVEGE and QVEGT respectively, Table 1) (Pistocchi et al. 2008 Scores from the BN model for water supply ranged potentially between -1 and +3 (Fig. 3), which were rescaled to the 0 -1 interval, with the new 0 placed at 0.25.Negative values (or values < 0.25 within the rescaled interval) would suggest a negative balance between water losses and water recharge: for instance, when losses within the soil and vegetation

LULC contribution to ES delivery
The (climate-based) results from the BN were then weighted according to the underlying land use and land cover (LULC) type, to better reflect the differences in the relative contribution of LULC types to support ES provision.We followed the 'ES potential matrix' method provided in Burkhard et al. (2014) which consists of plotting ES on the x-axis and geophysical spatial units (land cover types) on the y-axis, filled with integers from 0 to 5. Burkhard et al. (2014) distinguish between 'potentials' and 'flows', with the former indicating 'the potential maximum yield of selected ES' (Burkhard et al. 2012), and the latter indicating 'used set (bundles) of ES and other outputs from natural systems in a particular area within a given time period' (Burkhard et al. 2014).We chose to use the 'ES potential' since we do not have sufficient information to predict the actual demand for the selected ES into the future.The scores provided in Burkhard et al. (2014) are based on a typical central European landscape, classified according to the CORINE land cover classification system, and evaluated for the summer time, before harvest.To adapt the matrix to our study: • The rounded up average from the scores of Burkhard et al. (2014) was taken, when a LUISA LULC type included more than one CORINE land cover type; • All scores were rescaled to the 0.0 to 1.0 range.
Fig. 4 is a visual representation of the ES per LULC type scores, Suppl.material 2 shows the key between LUISA and CORINE.As a practical example, for an area of ca. 25 x 25 km (the land model resolution) with two LULC types (80% Arable land and 20% Forest), a score of 1 for 'Water flow regulation' resulting from the BN, would become 0.4 within the arable areas and 0.6 in forest areas (Fig. 4 and Suppl.material 2).The spatial resolution of the resulting map matched the LULC input, i.e. 1 ha (100 x 100 m grid cells).

Assessing relative changes in ES provision
To highlight areas where climate change is estimated to induce changes in the CLM inputs, we first extracted the percent change between 2050 and the baseline for each bio-physical parameter used within the BN model.We expected that major changes in the BN inputs would also propagate through the relevant ES.
Subsequently, following the ES per LULC matrix, we assessed the potential supply of ES, for the baseline and the 2050 scenario.From these two assessments we calculated the percent change over time and across geographical space, for each ES.The use of percentage also allowed us to compare changes in ES supply across ES see also Maes et al. 2015 ).
Further to this assessment, we extracted all areas where LULC was predicted to remain the same over time and we quantified the changes in ES provision, solely due to changing climatic conditions.This refinement allowed us to quantify the ES climate-induced changes in relation to those caused by changes in both drivers (LULC and climate).
Lastly, we highlighted the greatest positive and negative changes in ES provision occurring in forests throughout Europe.

Changes in the bio-physical variables supporting ecosystem processes
Climate change following the model predictions for the A1B scenario is expected to have a profound impact on the different ecological functions used in this study to derive ecosystem services (Fig. 5) The expected plant productivity, in general, will increase across Europe.This is visible from the expected changes in exposed one-side leaf area index (ELAI), and in the root biomass with increasing fine root carbon (FROOTC) and total soil organic matter carbon (TOTSOMC).The magnitude of these changes is relatively small, but widely spread across the study area.Aquifer recharge rates (QCHARGE) and the daily fire probability (FIRE_P) are expected to exhibit substantial changes in both directions (positive and negative), widely spread across the study area.Fire probability is expected to increase in north and south Europe while decreases are foreseen in the areas in between.Aquifer recharge rates are expected to fall in the Mediterranean, following the altered patterns of rainfall (RAIN).
Four functions, surface runoff (QCOVER), ground evaporation (QSOIL), canopy evaporation (QVEGE) and canopy transpiration (QVEGT) are expected to change in both directions for very limited regions, while for the majority of the regions small changes are foreseen.
Water in the unconfined aquifer (WA) is expected to undergo small and spatially localised changes while the canopy top (HTOP) will slightly change across Europe.

Changes in ES provision due to potential changes in LULC and climate
The next three assessments quantify the relative provision of ES for the baseline and for 2050, under a scenario of changes potentially occurring both drivers (LULC and climate).
Air quality regulation Fig. 6 shows air quality regulation resulting from the BN models over LULC types.Air quality regulation is displayed as a relative index between 0 and 1 for the baseline, and the 2050 and as a percent change between 2050 and the baseline.The index for air quality regulation had mean = 0.27 and SD = 0.24 for the baseline, and mean = 0.29 and SD = 0.26 for 2050.The percent change between 2050 and the baseline showed that most of the values were bound by the ±10% change with very few values towards the extreme ends (mean = 1.2 and SD = 11.6).The visual representation, therefore, was chosen to reflect this distribution.Areas of greatest negative change were mainly located in the EU-28 Scandinavian countries, in the Iberian Peninsula and on the South-East part of Europe.).These patterns were generally conserved over time (2050 vs. baseline), with the vast majority of the areas registering changes between -0.4 and + 2.5 % across the five regions (Fig. 7, panel 'Change over time').South Europe, however, also registered noticeable negative changes, over about 30% of the area, whilst Central and North Europe showed changes in scores ≥2.5% over ca.20% of their area.In addition, North Europe was the region with the greatest extent of the most negative changes.

Soil erosion control
Fig. 8 shows soil erosion control resulting from the BN models over LULC types.Soil erosion control is displayed as a relative index between 0 and 0.9 for the baseline and the 2050, and as a percent change between 2050 and the baseline.The index for erosion control had mean = 0.32 and SD = 0.26 for the baseline, and mean = 0.34 and SD = 0.27 for 2050.Global change impacts on ecosystem services: a spatially explicit assessment ...
The percent change between 2050 and the baseline showed that most of the values were within -2.8% and +6.8% change, with very few values towards the extreme ends (mean = 1.2 and SD = 13.4).The visual representation, therefore, was chosen to reflect this distribution.Extended areas of greatest negative change were located both in North and South Europe; in both regions, however, there were also noticeable areas of positive change.
Fig. 9 quantifies the distribution of erosion control scores and relative change over time, across the five European regions.North Europe was characterised by a rather homogeneous representation of all classes of scores, whilst in the BI the vast majority of the region was characterised by the lowest scores (≤0.2 over 60% of the region).In all other regions, ca.20% of the extent was characterised by high scores (EC ≥0.7).
The percent change over time (2050 vs. baseline) revealed that across all five regions the vast majority of their area was characterised by a little positive increase in the capacity of controlling soil erosion (although with similar SD, this effect is likely to be non-significant).Interestingly, North Europe was characterised by the greatest extent of both positive and negative changes (Fig. 9, panel 'Change over time').
Water flow regulation Fig. 10 shows water flow regulation resulting from the BN models over LULC types.Water flow regulation is displayed as a relative index between 0 and 0.7 for the baseline, between 0 and 0.8 for 2050, and as a percent change between 2050 and the baseline.The index for water flow regulation had mean = 0.25 and SD = 0.11 for the baseline, and mean = 0.24 and SD = 0.11 for 2050.The percent change between 2050 and the baseline showed that most of the values were within -3.4% and +2.4% change, with very few values towards the extreme ends (mean = -0.4 and SD = 6.3).The visual representation, therefore, was chosen to reflect this distribution.Extended areas of greatest negative change were located in the Northern-most part of the EU-28 Scandinavian countries as well as in the South-East part of Europe and    The percent change over time (2050 vs. baseline) showed that the largest extent of negative change was recorded in South Europe, immediately followed by North Europe.The eastern part of Central Europe was characterised by a rather homogenous change cross all six classes, whilst in the British Isles we observed the largest extent of positive change (Fig. 11, panel 'Change over time').In general, it is noticeable that the magnitude of change does not include extreme values: ad-hoc inspections revealed in fact that the greatest changes in absolute values were often caused by land conversion from one class to another one: for instance, if an area was predicted to undergo conversion from forest to arable land between the baseline and 2050, the score for the capacity of the underlying ecosystems to control soil erosion would change from 1 (forest) to 0.2 (arable land) (Fig. 4) (i.e.-80% change).

Changes in ES provision due to changes in climate alone
Air quality regulation Fig. 12 shows the change in air quality regulation over time (2050 vs. baseline), in areas where LULC was predicted to remain the same.The magnitude of change ranged between ±26.6% (mean = 0.7 and SD = 2.7).The regional breakdown shows that the largest extent of negative changes was predicted in the British Isles (BI).This is an interesting result which contrasts with the predictions obtained when including also changes in LULC (Fig. 7), where the results showed that the majority of the BI area would only undergo changes between -0.44 and + 2.5 %, whilst the most negative changes were observed in less than 5% of the BI area.Additionally, when looking only at areas with the same LULC over time, South Europe recorded the greatest extent of the most negative changes, while this was recorded in North Europe when also accounting for changes in LULC (Fig. 7).In all cases, therefore, looking only at areas that are predicted not to undergo LULC change can help parsing out the effects of climate from those of land conversion.

Soil erosion control
Fig. 13 shows the change in soil erosion control over time (2050 vs. baseline), in areas where LULC was predicted to remain the same.The magnitude of change ranged between +13.3% and -6.7% (mean = 0.7 and SD = 1.6).The regional breakdown shows that the largest extent of negative changes was predicted in the NE, where it was also recorded the largest extent of the most negative change.These results are consistent with the predictions that also accounted for LULC changes (Fig. 9) but they are greater in their extent; consequently, also the extent of positive changes in this region is substantially lower than the others.
Water flow regulation Fig. 14 shows the change in water flow regulation over time (2050 vs. baseline), in areas where LULC was predicted to remain the same.The magnitude of change ranged between +19.8% and -13.3% (mean = -0.6 and SD = 2.5).The regional breakdown shows that the largest extent of negative changes was predicted in South Europe, similarly to what already observed when accounting also for LULC.The largest extent of positive changes was recorded for the British Isles, again in accordance to what observed when also accounting for LULC.

Changes in forest areas
Fig. 15 shows the change in ES provision over time in forest areas, using quantile split computed from each of the ES percent change in forest areas only (i.e., masking forest from Figs 6, 8, 10).Quantiles were coded -2, -1, 0, +1, and +2.After summing all quantiles for the three ES, the results ranged from -6, indicating greatest decrease in ES provision over time, to +6, indicating greatest increase in ES provision over time.As we stated in the methods' description, only forest areas that were not predicted to change their extent over time were considered for this analysis: any change in ES provision, therefore, is the direct effect of change in climatic conditions affecting bio-physical processes underpinning the ES of interest.Since the scores attributed to the change are arbitrary, they have been translated into their more meaningful qualitative indicator featured in Fig. 15 (i.e., ranging from greatest negative change to greatest positive change).The map in Fig. 15 shows that, in general, most negative changes are expected in Southern European forests, and partially in the northernmost forests of Scandinavian countries.On the contrary, positive changes are expected in Central European forests, and in forests along the Atlantic coast of Spain and along the Pyrenees.

Discussion
With this study we tested the feasibility to link climate change and changes in land cover and land use to show changes in the supply of three regulating ecosystem services, across terrestrial ecosystems.Additionally, we identified areas where greatest positive and negative changes are expected, in a given ecosystem (here represented by a LULC type).

Main highlights
A few points can be highlighted from our work: • Mapping the predicted changes in the bio-physical functions underpinning the three regulating ES, helped us understand some of the predicted changes in the dependent ES.For instance, inspecting the changes in the components of water flow regulation (the ES with the predicted greatest negative change), showed noticeable changes for several of them: 'Aquifer recharge rate' (QCHARGE) above all, with great magnitude of change (positive and negative), widely spread across Europe; 'Water in the unconfined aquifer' (WA), with localised but important changes in the Southern part of Europe; and 'Precipitation' (RAIN), characterised by small magnitude of change, but widely spread across Europe.

•
Inspecting the changes in LULC, revealed potential major effects on ES supply.

•
We considered climate change and LULC as additive effects.But the changes in both drivers are likely to lead to more complex effects than expected by their additions (Elmhagen et al. 2015).Similarly, the combined increase of temperature, rainfall and CO is likely to lead to complex interactions and feedbacks: for instance, increased CO reduces stomatal conductance, which reduces transpiration and counteracts a potential increase in evapotranspiration caused by warming (Kløve et al. 2014).

The role of local management strategies
We performed the analysis at the extent of EU-28, using climate data with a spatial resolution of ca. 25 km and LULC information with a much finer spatial resolution (100 m).We did not include local management strategies which might be in place in certain areas, as this information cannot be readily integrated into our models and would require arbitrary decisions on conversion factors and weights, which could have a significant influence on the predictions.CC scenarios (IPCC 2000), soil erosion control and water provision along the Llobregat river basin (Spain) are going to be negatively affected; this is likely due to the increasing frequency of flood and drought events predicted across the Iberian Peninsula.
Refining the European-wide results with local scale information would allow us to account for processes that vary at very local scale, such as the effects of pressure on hydrology.
Local scale information would also allow us to evaluate the combined effects of different mitigation strategies: • Kiedrzyńska et al. ( 2015), for instance, showed how effective management of river floodplains can be used to mitigate against the predicted increase in flood events; they highlighted a combination of different mitigation strategies, involving modelling, land use management for better water retention capacity, and infrastructure.It is also widely recognised, however, that more information is needed on several mechanisms at the basis of ecosystem functioning (and hence, ES provision), such as groundwater recharge (Kløve et al. 2014).

Changes in ES supply in forests
When assessing the results from changes induced by climate-only, we chose forests as an ES provider as they support the provision of many ES, and certainly of all the ES assessed here (Fig. 4).
• Due to constraints on the input data, we only considered a broad definition of forest.
A differentiation between forest types (e.g.broad-leaves, coniferous, etc.) is not currently contemplated by LUISA.A workaround, which however was not feasible with the available resources, is to use the CORINE classification to differentiate between different types of forest, and projected climate change to find the closest analogous climate for a set of representative woody species.This approach, however, is poor for at least two reasons: (i) it assumes that the closest suitable climate will be within species' dispersal capacity, contradicting current literature (Ohlemüller et al. 2006a, Ohlemüller et al. 2006b) ; and (ii) it does not account for human interventions reacting to climate change or policy directives.

•
Forests are expected to increase under the climatic projections of the next decades, for instance due to increased availability of CO and, at least in Europe, as a consequence of agricultural land-abandonment and policy strategies.Yet, an increased in forest extent does not necessarily imply increased biodiversity (particularly for mono-specific forest plantations grown mainly for biomass production).In addition, evidence shows that several climate change projections are expected to produce unsuitable conditions for many species, and closest climate analogues beyond species' dispersal capacity (Ohlemüller et al. 2006a, Ohlemüller et al. 2006b).These two elements do not support the idea of stable forest communities anymore, implying the need of societal interventions to ensure forest's viability.Therefore, the increased in ES supply provided by a certain ecosystem, might only apply to a few ES and would lead to biased conclusion if used as the only criterion to assess the effects of CC and LULC changes.
Final remarks Some final remarks, but not of less importance, are on the assumptions and limitations of our method:

•
We have only considered one CC scenario (A1B) due to constraints on the input data.Our results, therefore, offer two snapshots in time (baseline and 2050 projections) of a very complex and wider reality.

•
We have captured some of the uncertainty related to CC projections by considering three different Regional Climate Models and by accounting for their differences while building the Bayesian Network, as well as when summarising their predictions.But the complexity of the subject would as well grant an examination of the uncertainty related to each of its components, from the input data to the processes and, where applicable, choice of weights, from the existing knowledge to the available tools.As noticed for instance by Grêt-Regamey et al. (2013), however, 2 such an analysis would require a separate trans-disciplinary process, to be really informative.

•
We have assumed the same scaling properties between the ES considered in our study: as such, a unit change in one ecosystem was weighted equally when deriving the quantile split.In absence of sufficient bio-physical evidence we lack the necessary information to support different choices.

•
We have assumed that any perturbation will have the same effect on the system: in other words, we have not included 'tipping points', i.e. critical thresholds at which a tiny perturbation can qualitatively alter the state or development of a system (Lenton et al. 2008).Incorporating this information remains a major research challenge, already highlighted, for instance, by Ciscar Martinez et al. (2014).Yet it is a key element to determine planetary boundaries, which are defined upstream the tipping point to buffer uncertainty, and should be used to guide human choices towards sustainable development (Steffen et al. 2015).

Conclusions
This work represents a preliminary spatially explicit assessment of the CC impacts on ES provision.As anticipated by Ciscar Martinez et al. (2014), climate change effects on ES provision is a new area of investigation, characterised, amongst others, by complex interactions between different components of the system (e.g.climate and LULC, climate and ecosystem processes) and by knowledge gaps.These are clearly major challenges to address, yet the scale of global change requires prompt actions to mitigate or adapt to the new conditions.Further research is needed not only to expand the analysis to the ES not considered here, but also to incorporate processes and scaling properties of the systems considered as they become available, as well as to account for spatial dependencies.
Despite these limitations, our work has also some highlights: This work is based on semi-automated procedures which can be routinely repeated as more (or better) information becomes available.This is particularly important given that this area of research is expanding, and can therefore provide evidencebased scientific and technical support to the EU policy cycle.
In addition to these methodological highlights, the results allow us to draw at least two general conclusions: • Positive and negative changes are spatially distributed, with a tendency for negative changes to be distributed at the periphery of the study area (i.e.mainly in the northern and southern areas), and positive changes to be distributed towards the central parts.This is particularly noticeable when looking at the changes affecting water flow regulation.This general pattern is probably a consequence of climatic conditions being already at the extreme of their reference range in the northern and southern areas.• LULC change is a stronger driver than climate change on the provision of ES: this has important implications as it suggests that the impacts of climate change can be mitigated through adaptive ecosystem management.

Figure 2 .
Figure 2. Main data, processes and outputs to estimate changes in ES provision.When relevant, sources of data are indicated within round brackets () and defined within the main text.The spatial resolution of information from LUISA and CLM is specified within square brackets [].BN = Bayesian Network; ES = Ecosystem Service(s); LULC = Land use / land cover

2 2Figure 3 .
Figure 3. Schematic representation of the BN built for the assessment of the selected ES (Soil erosion control, Air quality regulation, Water flow regulation).Bio-physical components from the CLM are indicated within square brackets [].See Table1for their definition and see the main text for a description of each relevant ES.

Figure 4 .
Figure 4. Relative capacity (y-axis) of LULC types (x-axis) to supplying regulating ES (air quality regulation, soil erosion control and water flow regulation).LULC types are defined as in LUISA (Baranzelli et al. 2014), scores are based on a modified version of the ecosystem service potential matrix of Burkhard et al. (2014).See main text for additional details.Abnd: Abandoned; Trnstnl: Transitional.
. Positive changes are shown in red shades, whilst negative changes are shown in blue shades; no change (or very little change in relation to the relevant magnitude) is displayed in yellow.The changes are depicted for the period between 2050 and the baseline.The magnitude and the spatial extent of change vary greatly across different variables.

Figure 5 .
Figure 5. CLM model outcomes showing climate-induced percent change between 2050 and the baseline, for the bio-physical variables used as inputs of the ES-BN models.See Table 1 for a definition of the variable names.The magnitude and spatial extent of change can be categorised in five main patterns: (1) Great magnitude of change, widely spread across the study area (FIRE_P and QCHARGE); (2) Great magnitude of change for very limited regions, and the majority of the regions laying around small changes (QOVER, QSOIL, QVEGE and QVEGT); (3) Small magnitude of change, but widely spread across the study area (e.g.ELAI, FROOTC, RAIN, TOTSOMC); (4) Small magnitude of change, spatially localised (WA); (5) Very small change, but widely spread (HTOP).Distribution maps are displayed using nonprojected Geographic Reference System (WGS84 Datum).

Figure 6 .
Figure 6.Air quality regulation for baseline, 2050 and as a percent change over time (2050 vs. baseline).For the baseline and the 2050 maps, classes represent equal intervals over the range of rendered predictions; for the percent change over time, the visual representation reflects the distribution of the values.Maps are displayed using the Projected Reference System LAEA.See main text for additional details.

Fig. 7
Fig.7quantifies the distribution of air quality regulation scores and relative change over time, across five European regions (Suppl.material 3) identified on the basis of their climatic characteristics(Ciscar Martinez et al. 2011).In general, scores ≤0.2 characterised at least 40% of each region's extent (a little less for North Europe, NE).This percentage was much greater in South Europe (SE) and on the British Isles (BI), where values of AQ ≤0.2 extended up to and above 60% of each region, respectively.The North part of Central Europe (CEN), recorded the greatest percentage in high supply of air quality regulation (scores ≥0.8).These patterns were generally conserved over time (2050 vs. baseline), with the vast majority of the areas registering changes between -0.4 and + 2.5 % across the five regions (Fig.7, panel 'Change over time').South Europe, however, also registered noticeable negative changes, over about 30% of the area, whilst Central and North Europe showed changes in scores ≥2.5% over ca.20% of their area.In addition, North Europe was the region with the greatest extent of the most negative changes.

Figure 7 .
Figure 7. From left to right, air quality regulation for the baseline, the 2050 and as percent change over time.The scores are shown in relation to the percent extent of five European regions, defined after grouping the 28 Member States (MS) according to their climatic characteristics (Ciscar Martinez et al. 2011).The chart on the right shows the percent change over time in the predicted air quality regulation, for the same set of regions.In all charts, the y-axis shows the different regions (SE: South Europe; CEN: Central Europe North; CES: Central Europe South; NE: North Europe; BI: British Isles).See Suppl.material 3 for a grouping of MS within each region.See main text for additional details.

Figure 8 .
Figure 8. Soil erosion control for the baseline, the 2050 and as a percent change over time (2050 vs. baseline).For the baseline and the 2050 maps, classes represent equal intervals over the range of rendered predictions; for the percent change over time, the visual representation reflects the distribution of the values.Maps are displayed using the Projected Reference System LAEA.See main text for additional details.

From
left to right, soil erosion control for the baseline, the 2050 and as percent change over time.The scores are shown in relation to the percent extent of five European regions, defined after grouping the 28 Member States (MS) according to their climatic characteristics (Ciscar Martinez et al. 2011).The chart on the right shows the percent change over time in the predicted soil erosion control, for the same set of regions.In all charts, the y-axis shows the different regions (SE: South Europe; CEN: Central Europe North; CES: Central Europe South; NE: North Europe; BI: British Isles).See Suppl.material 3 for a grouping of MS within each region.See main text for additional details.

Figure 10 .
Figure 10.Water flow regulation for the baseline, the 2050 and as a percent change over time (2050 vs. baseline).For the baseline and the 2050 maps, classes represent equal intervals over the range of rendered predictions; for the percent change over time, the visual representation reflects the distribution of the values.Maps are displayed using the Projected Reference System LAEA.See main text for additional details.

Fig. 11
Fig.11quantifies the distribution of water flow regulation scores and relative change over time, across the five European regions.For the baseline, the vast majority of the area within each region was characterised by scores between 0.2 and 0.3, with the exception of North Europe, where the vast majority of the region's extent scored 0.4.The British Isles and North Europe also recorded noticeable areas within the interval 0.5 -0.6.The 2050 predictions showed, again for North Europe, the vast majority of the region's extent within the interval 0.4 -0.5; whilst areas characterised by the second highest score (0.6) were mainly confined to the British Isles.

Figs 12 ,
Figs 12, 13, 14  show the percent change in ES provision for air quality regulation, soil erosion control, and water flow regulation respectively, in areas where the only driver of change was climate (in other words, where land use land cover (LULC) was predicted to remain unchanged over time).

Figure 11 .
Figure 11.From left to right, water flow regulation for the baseline, the 2050 and as percent change over time.The scores are shown in relation to the percent extent of five European regions, defined after grouping the 28 Member States (MS) according to their climatic characteristics (Ciscar Martinez et al. 2011).The chart on the right shows the percent change over time in the predicted air quality regulation, for the same set of regions.In all charts, the y-axis shows the different regions (SE: South Europe; CEN: Central Europe North; CES: Central Europe South; NE: North Europe; BI: British Isles).See Suppl.material 3 for a grouping of MS within each region.See main text for additional details.

Figure 12 .
Figure 12.Change in air quality regulation over time (2050 vs. baseline), in areas where LULC was predicted to remain unchanged.Areas mapped in dark red indicate largest percent decrease, areas mapped in dark green indicate the largest percent increase.The bar chart shows the change as a relative extent percent, for five European regions defined after grouping the 28 Member States according to their climatic characteristics (Ciscar Martinez et al. 2011).The map is displayed using the Projected Reference System LAEA.

Figure 13 .
Figure 13.Change in soil erosion control over time (2050 vs. baseline), in areas where LULC was predicted to remain the same.Areas mapped in dark red indicate largest percent decrease, areas mapped in dark green indicate the largest percent increase.The bar chart shows the change as a relative extent percent, for five European regions defined after grouping the 28 Member States according to their climatic characteristics (Ciscar Martinez et al. 2011).The map is displayed using the Projected Reference System LAEA.

Figure 14 .
Figure 14.Change in water flow regulation over time (2050 vs. baseline), in areas where LULC was predicted to remain the same.Areas mapped in dark red indicate largest percent decrease, areas mapped in dark green indicate the largest percent increase.The bar chart shows the change as a relative extent percent, for five European regions defined after grouping the 28 Member States according to their climatic characteristics (Ciscar Martinez et al. 2011).The map is displayed using the Projected Reference System LAEA.

Figure 15 .
Figure 15.Change in ecosystem service provision over time (2050 vs. baseline) in forest areas over EU-28 (left) and close-up examples (right).Dark shades of green correspond to great positive change, and indicate areas with the greatest increase in ES provision, dark shades of red correspond to great negative change and indicate areas with the greatest decrease in ES provision.Maps are displayed using the Projected Reference System LAEA.

Table 1 .
Bio -physical variables from the CLM (proxies for ecosystem functions, EF, or ecosystem processes, EP) and related ecosystem services (ES).Positive relationships between EF / EP and ES are indicated with '+', negative relationships with '-'.Global change impacts on ecosystem services: a spatially explicit assessment ...

[
Eq. 3]layers (ground evaporation and evapotranspiration) exceed the recharge provided by rainfall and/or the water reserve available within the water table.See section 'LULC contribution to ES delivery' for further details on the relative contribution of LULC classes to water flow regulation.Equation 3 summarises the model: Our outcome, therefore, offers an overview of what can be expected across EU-28.Hence, patterns emerging from our work should be compared with, and complemented by, the increasing number of studies that are becoming available both at European and local scale: • At European scale for instance, Panagos et al. (2015b) recently proposed a modified version of the RUSLE model (RUSLE2015), where cover-management factor (Panagos et al. 2015a) and support practices in arable land are incorporated into their model, to estimate soil loss rates by water for the reference year 2010.Despite our approaches and input variables are different, the predicted outcomes show some similarities: for instance, general low rates of soil erosion throughout

•
To our knowledge, this is the first spatially explicit example of an EU-wide assessment linking CC, LULC and ES.It brings together scientifically accepted models (e.g.CLM and LUISA) and the latest tools for ES-LULC mapping (e.g. the ES potential matrix), to derive the potential changes in ES provisions, caused by CC alone, as well as by the combined change in climate and LULC.• Linking climate and LULC models enabled us to produce results at a spatial resolution finer than what would be obtained using climate models only.This is important not only to better represent variation and processes acting at local scale, but also to better inform local scale interventions.•