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 climate-induced 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.

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 non-overlapping 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:

• Air quality regulation
• Soil erosion control
• 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.

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).
• 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:

Air quality regulation = Rain (w) + Vegetation control – Fire probability [Eq. 1]

Where:

• Rain (w) = 0.8 * RAIN;
• Vegetation control = Foliage (ELAI) + Canopy height (HTOP)

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 2 summarises the model:

Erosion Control = Protection – Surface runoff [Eq. 2]

Where:

• Protection = Vegetation (FROOTC)+ Soil organic matter (TOTSOMC)
• Surface runoff = QOVER

​ 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, Vigerstol and Aukema 2011).
• Ground water: We used water in the unconfined aquifer and aquifer recharge rate (WA and QCHARGE respectively, Table 1) to represent groundwater (Jujnovsky et al. 2012).

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 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:

Water supply = Water in the unconfined aquifer + Aquifer recharge rate + Available water [Eq. 3]

Where:

• Water in the unconfined aquifer = WA
• Aquifer recharge rate = QCHARGE
• Available water = RAIN – Water loss;
• Water loss = Ground evaporation (QSOIL) + Canopy evaporation (QVEGE)+ Canopy transpiration (QVEGT)

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).

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.

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). 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.

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.

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.

Fig. 7 quantifies 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.

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.

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 map for the baseline water flow regulation did not show a spatial pattern in the distribution of scores, whilst the 2050 prediction showed some regions with greater capacity of regulating water flow (e.g. Northern part of the British Isles, vast areas in the EU-28 Scandinavian countries, Northern part of The Netherland and of Portugal, to give a few examples).

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 Iberian Peninsula. Areas of positive change were visible in the Southern part of Central Europe, in the Northern part of the British Isles and of Portugal.

Fig. 11 quantifies 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.

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’).

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).

Compared to Figs 6, 8, 10 (ES change percent over time), these new maps only include areas where changes in ES provision are caused exclusively by climate-driven changes in the underlying ecosystem bio-physical processes.

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).

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.