Climate change, reforestation/afforestation, and urbanization impacts on evapotranspiration and streamflow in Europe

Since the 1950s, Europe has undergone large shifts in climate and land cover. Previous assessments of past and future changes in evapotranspiration or streamflow have either focussed on land use/cover or climate contributions or on individual catchments under specific climate conditions, but not on all aspects at larger scales. Here, we aim to understand how decadal changes in climate (e.g. precipitation, temperature) and land use (e.g. deforestation/afforestation, urbanization) have impacted the amount and distribution of water resource availability (both evapotranspiration and streamflow) across Europe since the 1950s. To this end, we simulate the distribution of average evapotranspiration and streamflow at high resolution (1 km²) by combining (a) a steady-state Budyko model for water balance partitioning constrained by long-term (lysimeter) observations across different land use types, (b) a novel decadal high-resolution historical land use reconstruction, and (c) gridded observations of key meteorological variables. The continental-scale patterns in the simulations agree well with coarser-scale observation-based estimates of evapotranspiration and also with observed changes in streamflow from small basins across Europe. We find that strong shifts in the continental-scale patterns of evapotranspiration and streamflow have occurred between the period around 1960 and 2010.

In much of central-western Europe, our results show an increase in evapotranspiration of the order of 5 %–15 % between 1955–1965 and 2005–2015, whereas much of the Scandinavian peninsula shows increases exceeding 15 %. The Iberian Peninsula and other parts of the Mediterranean show a decrease of the order of 5 %–15 %. A similar north–south gradient was found for changes in streamflow, although changes in central-western Europe were generally small. Strong decreases and increases exceeding 45 % were found in parts of the Iberian and Scandinavian peninsulas, respectively. In Sweden, for example, increased precipitation is a larger driver than large-scale reforestation and afforestation, leading to increases in both streamflow and evapotranspiration. In most of the Mediterranean, decreased precipitation combines with increased forest cover and potential evapotranspiration to reduce streamflow. In spite of considerable local- and regional-scale complexity, the response of net actual evapotranspiration to changes in land use, precipitation, and potential evaporation is remarkably uniform across Europe, increasing by ∼ 35–60 km³ yr⁻¹, equivalent to the discharge of a large river. For streamflow, effects of changes in precipitation (∼ 95 km³ yr⁻¹) dominate land use and potential evapotranspiration contributions (∼ 45–60 km³ yr⁻¹). Locally, increased forest cover, forest stand age, and urbanization have led to significant decreases and increases in available streamflow, even in catchments that are considered to be near-natural.

1 Introduction

Streamflow provides an integrated signal in both space and time of all upstream changes in the terrestrial hydrological cycle. At smaller timescales of days to weeks, streamflow reflects the weather conditions and precipitation in the recent past. At longer (multi-year) timescales, over which internal catchment storage changes become much smaller than the amount of water passing through the catchment system, streamflow reflects the amount of water that passes through aquifers and dams (the “water yield”), which is the portion of precipitation that is not returned to the atmosphere via evapotranspiration. The water yield represents the average water flux that can potentially be exploited for human benefit in a sustainable way. Quantifying and understanding past and future changes in water availability in rivers and groundwater systems, reflecting the integrated signal of all net changes in the water cycle upstream, are not only of key importance to water resource management and planning, but are also a major scientific challenge given the uncertainties and limitations in both observations and models (Wang, 2014). This is in particular true for Europe, where strong changes in land use (in particular urbanization, reforestation and afforestation; see Fuchs et al., 2013) and climate (van der Schrier et al., 2013; Caloiero et al., 2018; Bach et al., 2018) have occurred since the 1960s.

Several studies have focussed on large-scale changes in evapotranspiration and/or water availability. In one of the first global studies, Milly et al. (2005) analysed climate-driven changes in water availability from an ensemble of climate models and found a general drying of transitional regions and a wetting of current humid and colder regions. Over Europe, the study reported a strong latitudinal gradient in average water fluxes increasing in strength from the 20th to the 21st centuries, with decreasing availability trends in the Mediterranean and increasing trends in northern Europe. Gerten et al. (2008) showed that globally, precipitation changes were the biggest drivers of changes in runoff, but land use change also had a considerable effect. Changes in Northern Hemisphere streamflow over the past decades have likely also been impacted by decadal changes in solar radiation (the so-called global “dimming” and “brightening”; see Teuling et al., 2009; Gedney et al., 2014). Other studies have focussed on the impact of anthropogenic land cover change on evapotranspiration. Sterling et al. (2012) found a 5 % reduction in global evapotranspiration due to land cover conversion, resulting in a 7.6 % increase in global average streamflow. Other studies have highlighted strong decadal-scale variability in global average evapotranspiration over the recent decades related to El Niño–Southern Oscillation (Jung et al., 2010; Miralles et al., 2014). In spite of the direct link between average evapotranspiration and streamflow, few studies have addressed changes in both fluxes simultaneously.

Because streamflow is impacted by many factors, which often have opposing effects, changes in streamflow should be considered at small scales at which individual factors can be understood and quantified rather than at larger river basin scales. Although several long discharge records exist for large river basins, changes that occur at the sub-basin level are often obscured by opposing effects in other parts of the basin. In a landmark study, Stahl et al. (2010) addressed this limitation by analysing streamflow changes in Europe from a dataset of relatively small river basins with limited human influence. They reported a diverging pattern of streamflow trends over the past decades, with negative trends in annual mean streamflow in many parts of the Mediterranean and central Europe and predominantly positive trends in western Europe and parts of Scandinavia. While the longer-term and long-range variability of streamflow in these basins and its relation to circulation indices is generally well understood at the interannual and decadal timescales (Gudmundsson et al., 2011; Hannaford et al., 2013), significant uncertainty exists in understanding the regional-scale variability in trends since these are not well reproduced by the current generation of hydrological models (Stahl et al., 2012). Previous regional case studies across Europe have reported a sensitivity of long-term water balance partitioning to both climate and land use change (Parkin et al., 1996; van Roosmalen et al., 2009; van der Velde et al., 2013; Renner et al., 2014; Pijl et al., 2018). Thus, quantifying changes in streamflow requires accounting for changes in climate (precipitation and potential evapotranspiration) as well as changes in land use and/or land cover (Stonestrom et al., 2009). But whereas assessing the impact of climate on average streamflow is relatively straightforward, the role of land cover requires a more careful consideration.

At the smaller scale, land use, in particular forest cover, has long since been known to have a strong impact on average streamflow or water yield, with forested catchments having a much lower water yield compared to non-forested catchments (Bosch and Hewlett, 1982; Zhang et al., 2001; Brown et al., 2005; Farley et al., 2005; van Dijk and Keenan, 2007; Filoso et al., 2017). Based on analysis of paired catchment observations, a large majority of studies have found that removal of forest leads to an increase in water yield. While this is likely linked to higher average evapotranspiration over forest, the reverse has been reported for dry and warm summer conditions based on eddy-covariance observations from FLUXNET (Teuling et al., 2010). Somewhat surprisingly, average evapotranspiration rates for forested FLUXNET sites are generally slightly lower than for non-forested sites (Williams et al., 2012), which is seemingly inconsistent with other studies (e.g. Zhang et al., 2001), where annual evapotranspiration was inferred from the water balance (the so-called “forest evapotranspiration paradox”; see Teuling, 2018). A possible explanation for this discrepancy is the role of interception (van Dijk et al., 2015). Several studies (e.g. Gash et al., 1980; Zimmermann et al., 1999) have shown that interception can constitute a major term in the water balance of forested ecosystems, in particular in humid conditions (Calder, 1976; Ramírez et al., 2018). Controlled experiments on large non-weighable lysimeters covered with forest have shown that growing forest strongly reduces the water yield (Tollenaar and Ryckborst, 1975; Harsch et al., 2009; Müller, 2009; Teuling, 2018) and that this effect is somewhat larger for coniferous than for deciduous species. This is in line with results from a large number of basins in Sweden, where increases in forest cover and biomass (age) were the main factors explaining observed trends in inferred evapotranspiration (Jaramillo et al., 2018). This shows that forest cover area but also stand age need to be taken into account when evaluating land use change effects on water balance partitioning.

In contrast to forest cover, few studies have quantified the effects of urban area and urbanization on the long-term water balance partitioning. Runoff from urban areas is typically measured with a focus on short-term dynamics and event runoff ratios (Berthier et al., 1999) or runoff produced by impervious areas only (Boyd et al., 1993; Shuster et al., 2005). Evapotranspiration from urban areas, on the other hand, is typically measured or analysed over individual elements that make up the urban landscape, such as (un)paved areas (Ramamurthy and Bou-Zeid, 2014), green roofs, or trees (Pataki et al., 2011). Few studies have measured evapotranspiration at the urban landscape scale. In a study comparing measurements made over the Dutch cities of Rotterdam and Arnhem, Jacobs et al. (2015) found evapotranspiration rates to be generally low and to quickly drop in the days following rainfall, reflecting a strongly water-limited system. Similar results were found for the Swiss city of Basel (Christen and Vogt, 2004). This suggests that urban areas, because of their limited capacity to store water, might have much lower evapotranspiration and as a result might generate much higher streamflow than other land use types. This was also reported by DeWalle et al. (2000) based on statistical analysis of the long-term streamflow record in the United States. They found strong increases in streamflow in areas with heavy urbanization, which was attributed to a decrease in evapotranspiration.

In order to isolate and/or attribute the hydrological impact of climate change from that of changes in land use, different methods exist (see reviews by Wang, 2014; Dey and Mishra, 2017). The methods can be categorized into experimental approaches, hydrological modelling, conceptual approaches, and analytical approaches (Dey and Mishra, 2017). Typically, hydrological or land surface models run at hourly or daily resolution are used (Bosmans et al., 2017; Breuer et al., 2009; Viney et al., 2009; Dwarakish and Ganasri, 2015; Pijl et al., 2018). Such models often contain a high number of poorly constrained parameters and parameterizations, leading to large uncertainty in trend estimates (Arnell, 2011) or even disagreement in the direction of simulated trends (Melsen et al., 2018). When the research focus is on robust simulation of long-term rather than short-term changes, low-dimensional models with well-constrained parameters often perform well (Choudhury, 1999; Zhang et al., 2008). The Budyko model (Budyko, 1974) is an example of such a conceptual approach which allows for evaluation of combined land use and climate impacts on water availability (see, for example, Jiang et al., 2015). In spite of its extreme simplicity (parameterizations typically have only one parameter reflecting land surface characteristics), it has been applied successfully in numerous studies focussing on different controls on long-term water balance partitioning (Zhang et al., 2004, 2016; Roderick and Farquhar, 2011; Xu et al., 2013, 2014; Greve et al., 2014; Creed et al., 2014; Jiang et al., 2015; Wei et al., 2018). Although it is generally applied at coarse global grid resolution or to large river basins, other studies (e.g. Zhang et al., 2004; Redhead et al., 2016) have found the model to also work well for smaller basins or grid cells (<10 km²). This opens up the possibility for robust and parsimonious modelling of hydrological impacts at high spatial resolution.

The strong impact of land use on water balance partitioning at smaller scales, combined with the large-scale land use changes that have occurred over Europe over the past decades, leads to the question how they have impacted changing patterns of evapotranspiration and streamflow. Previous assessments of past and future changes in water balance partitioning have either focussed on land use/cover (Sterling et al., 2012) or climate contributions (Wilby, 2006; Gardner, 2009; Hannaford et al., 2013) or have focussed on smaller catchments under particular climate conditions (van Roosmalen et al., 2009; Renner et al., 2014; Pijl et al., 2018). Therefore, we aim to understand how recent decadal changes in climate (e.g. precipitation, temperature) and land use (deforestation/afforestation, urbanization) have impacted the amount and distribution of water resource availability across Europe since the 1950s. We address the hypothesis that land cover changes play a much more important role at the European scale than previously reported, even in basins which are assumed to have a limited human influence on the water cycle. To this end, we simulate the distribution of evapotranspiration and streamflow at high resolution (1 km²) by combining (a) a steady-state Budyko model for water balance partitioning constrained by long-term observations across different land use types, (b) a novel decadal high-resolution historical land use reconstruction, and (c) gridded observations of key meteorological variables. Simulations will be evaluated against state-of-the-art observation-based assessments of evapotranspiration and observed changes in streamflow.

2 Methods and data

Central to our approach is the formulation of the Budyko model as used by Zhang et al. (2004). As with any Budyko approach, it follows the central assumptions that the fraction of precipitation that returns to the atmosphere as evapotranspiration ET depends on the ratio between the average potential evapotranspiration PET and average precipitation P, rather than on their absolute values, and that a catchment's ET, when a catchment is subjected to a range of climate conditions, follows a single path in the ET ∕ P, PET ∕ P space. Good fits with observations at several spatial scales show that this assumption is generally justified. In the work by Zhang et al. (2004), the following equation was proposed for the dependency of ET∕P on PET∕P:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{ET}}{P}}=\mathrm{1}+{\displaystyle \frac{\mathrm{PET}}{P}}-{\left[\mathrm{1}+{\left({\displaystyle \frac{\mathrm{PET}}{P}}\right)}^w\right]}^{\mathrm{1}/w},\end{array}$$

in which w is a model parameter which is typically linked to catchment and/or vegetation properties (Li et al., 2013). Zhang et al. (2004) found w=2.63 to best fit observations for Australian catchments, with slightly lower values for grassed (w=2.55) and higher for forested catchments (w=2.84). While these different values confirm that w depends on land surface characteristics, the magnitude of this dependency at the scale of individual land use elements, rather than catchments with a land use mixtures of varying degrees, is probably larger. Based on analysis of remotely sensed Normalized Difference Vegetation Index (NDVI) and gridded global fields of ET, PET, and P at the 0.5^∘ resolution, Greve et al. (2014) reported values of 3.05 for grid cells with an NDVI of around 0.8, whereas grid cells with an NDVI of around 0.2 were found to follow w=1.63. In a similar study but using observed streamflow rather than estimated ET, Li et al. (2013) found w to depend on the basin-average fractional vegetation cover M according to $w=\mathrm{2.36}\times M+\mathrm{1.16}$. These studies show that w can show considerable variation even at relatively coarse scales.

In order to get the most realistic values for w for application at smaller scales (∼1 km²) at which land use is often fairly homogeneous and the effects on water balance partitioning are most pronounced, we constrain w by the best available observations for different land use types and made under European climate conditions. It should be noted that widely available FLUXNET observations are not used in this study, because they might show the opposite land use ET signal from water balance-based observations (the so-called forest evapotranspiration paradox; see Teuling, 2018). The latter are assumed here to be more reliable for our application. The observations used in this study primarily come from the long-term lysimeter stations, such as the ones at Rietholzbach (Seneviratne et al., 2012), St. Arnold (Harsch et al., 2009), Brandis (Haferkorn and Knappe, 2002), Eberswalde–Britz (Müller, 2009), Castricum (Tollenaar and Ryckborst, 1975), and Rheindahlen (Xu and Chen, 2005), several of which were also analysed in a previous study by Teuling (2018). These data are complemented by observations from a natural lysimeter at Plynlimon (Calder, 1976) under more humid climate conditions and flux observations made over the cities of Basel (Christen and Vogt, 2004), Arnhem, and Rotterdam (Jacobs et al., 2015). Long-term data are preferred to minimize impacts of interannual storage variations (Istanbulluoglu et al., 2012). By relying on lysimeter observations to constrain our Budyko parameters, we implicitly assume lysimeters (area varies from 1 to 625 m² for the larger lysimeters at Castricum) to behave similarly to landscape elements of 10⁶ m² (our grid cell size). The data are shown in Fig. 1 and listed in Table 1. It should be noted that the stations are not distributed evenly across Europe, but are mainly constrained to central-western Europe (Fig. A1).

Figure 1 Climate and land use controls on water balance partitioning from long-term flux observations. See Table 1 for origin of data points. The error bar indicates the total spread over multiple lysimeters at the Brandis site with different soil types (Haferkorn and Knappe, 2002). Curves are based on Eq. (3). Note that symbol shape indicates land cover, but colours indicate stand age in the case of forest. The dashed grey line indicates the energy limit for non-adjusted PET.

Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Jacobs et al. (2015)Christen and Vogt (2004)Jacobs et al. (2015)Harsch et al. (2009)Harsch et al. (2009)Harsch et al. (2009)Haferkorn and Knappe (2002)Xu and Chen (2005)Seneviratne et al. (2012)Seneviratne et al. (2012)Seneviratne et al. (2012)Seneviratne et al. (2012)Harsch et al. (2009)Harsch et al. (2009)Harsch et al. (2009)Harsch et al. (2009)Harsch et al. (2009)Harsch et al. (2009)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Tollenaar and Ryckborst (1975)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Müller (2009)Calder (1976)

Table 1 Data used in the Budyko analysis. Units of fluxes are in mm yr⁻¹.

¹Derived from E-OBS (P) and CRU (PET). ² Mean of 24 lysimeters listed, minimum value 478 mm yr⁻¹, and maximum 614 mm yr⁻¹ also shown as an error bar in Fig. 1. ³ Values digitized from Calder (1976).

Due to the smaller scale than applied in previous Budyko analyses, we initially find many points, in particular observations from forested lysimeters, to be located above the energy limit (grey dashed line in Fig. 1). This indicates that the long-term average yearly evapotranspiration (ET) exceeds the average potential evapotranspiration (PET). This is possible, for instance, due to evaporation of interception water by energy not captured in the formulation of PET (van Dijk et al., 2015). Therefore, we correct for the underestimation by introducing a so-called adjusted potential evapotranspiration (aPET) which is assumed to be proportional to the potential evapotranspiration and accounts for all processes affecting yearly ET for tall vegetation (including evaporation of intercepted water through advection) so that ET generally will not exceed aPET even for forested areas:

$$\begin{array}{}\text{(2)}& \mathrm{aPET}=c\times \mathrm{PET},\end{array}$$

resulting in the following expression for the Budyko curve:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{ET}}{P}}=\mathrm{1}+{\displaystyle \frac{\mathrm{aPET}}{P}}-{\left[\mathrm{1}+{\left({\displaystyle \frac{\mathrm{aPET}}{P}}\right)}^{w^{*}}\right]}^{\mathrm{1}/w^{*}},\end{array}$$

in which w^* is the value for w when aPET rather than PET is used. aPET should thus be interpreted as the maximum total yearly evapotranspiration that would occur under given climate conditions (PET ∕ P) and under optimal vegetation conditions (i.e. vegetation that is most efficient in returning precipitation to the atmosphere, in this case needleleaf forest), rather than a land use-specific crop factor. It should be stressed that this scaling is only done to move most observations within the energy and water limits in the Budyko space (so that we can obtain a fit with Eq. 1) and that it has no other impact on the results since the model is subsequently forced with aPET rather than PET. The resulting values for w^* that match the (lysimeter) observations are shown in Fig. 1. It was found that c=1.6 was required to ensure most observations would be located on the right-hand side of the energy limit ($\mathrm{ET}/P=\mathrm{PET}/P$). Subsequently, in all analysis we replaced PET with aPET, including Eq. (1), but also in the atmospheric forcing fields. It should be noted that while this procedure results in lower values for w^* that cannot be directly compared to values for w reported in previous studies, most of the simulated ET values are identical to the ones that would be simulated with the original model. We find the highest w^* for full-grown forest, indicating that any change towards this state due to reforestation or afforestation will increase ET given the same climate (P and PET). We distinguish between young stands (age <10 years, assumed to behave similarly to croplands/grasslands based on the data in Fig. 1 with $w^{*}=\mathrm{1.7}$), intermediate (age 10–20 years, $w^{*}=\mathrm{2.3}$), and older stands (age >20 years, $w^{*}=\mathrm{3.1}$); see also Fig. 1. In this way, we implicitly account for effects of increasing biomass, tree height, and stand age on ET and water yield reported in previous studies (Harsch et al., 2009; Jaramillo et al., 2018; Teuling, 2018). Conversely, urban areas have a low w^* of 1.3, indicating that urbanization will generally decrease ET. Finally, the long-term average streamflow or water yield at the pixel level is calculated from the catchment water balance,

$$\begin{array}{}\text{(4)}& Q\approx P-\mathrm{ET},\end{array}$$

under the assumption that storage changes (such as snow, soil moisture, or groundwater) and net lateral groundwater inflow/outflow can be neglected at the decadal (10-year) timescale. This timescale is chosen to align with the temporal resolution of the land use dataset and to minimize possible impacts of storage changes.

As input to our model as described above (Eqs. 2–4), we use gridded datasets of land cover and meteorological observations. All calculations were performed at a 1×1 km spatial resolution, which were later rescaled to a coarser resolution for visualization purposes. Historic land change information is based on the HIstoric Land Dynamics Assessment (HILDA, v2.0) model reconstruction of historic land cover/use change (Fuchs et al., 2013, 2015a, b). This data-driven reconstruction approach used multiple harmonized and consistent data streams such as remote sensing, national inventories, aerial photographs, statistics, old encyclopedias, and historic maps to reconstruct historic land cover at a 1×1 km spatial resolution for the period 1900 to 2010 in decadal time steps. The reconstruction provides information for six different land cover/use categories: forest, grassland (including pastures, natural grasslands, and shrublands), cropland, settlements/urban, water bodies and other (bare rock, glaciers, etc.). Here we only use the forest, grassland/cropland, and settlement classes. The reconstruction considers gross land changes, the sum of all area gains and losses that occur within an area and time period, unlike other reconstructions that focus on net changes only, calculated by area gain minus the area losses. Details on the net versus gross changes can be found in Fuchs et al. (2015a). The gross changes are used to derive forest stand age. Previous research has shown that not accounting for gross land use changes in reconstruction led to serious underestimations in the amount of total land use changes that have occurred (Fuchs et al., 2015a). The E-OBS v18 gridded dataset (Haylock et al., 2008) of observed precipitation at 0.25^∘ resolution and the CRU TS v4.02 gridded dataset (Harris et al., 2014) of observed potential evapotranspiration at 0.5^∘ resolution were used to force the model (Eqs. 2–4). Based on the joint availability of the HILDA, CRU, and E-OBS datasets, we selected two 10-year periods which were considered for analysis: 1955–1965 and 2005–2015. In the following, we will refer to these periods as 1960 and 2010 for simplicity. While the 10-year periods are often considered short for climate change detection, they resulted from a need to balance robust estimation of the mean climate without averaging out much of the underlying changes in both climate and land use. Changes over the intermediate 10-year periods were analysed, but since the trends were found to be mostly monotonic the results are not shown here (except for validation purposes in Fig. 5).

Figure 2 Climate characteristics over the period 1960–2010. Left panels show the mean precipitation and potential evapotranspiration (a and c, respectively), while the right panels indicate the change over the period 1960–2010 for precipitation (b) and potential evapotranspiration (d).

Figure 3 Land cover characteristics over the period 1960–2010. Left panels show the mean forest cover and urban fraction in 1960 (a and c, respectively), while the right panels indicate the change between the periods 1960 and 2010 for forest cover (b) and urban area (d).

Model simulations are validated and compared against observed yearly average streamflow changes in near-natural catchments and observation-based average evapotranspiration. The relative streamflow changes for the period 1962–2004 (normalized by the standard deviation of yearly streamflow) were used as presented in Stahl et al. (2010, their Fig. 2). Average evapotranspiration was derived from GLEAM v3.2a (Martens et al., 2017). The contribution of P, PET, and land use (through w^*) was assessed by performing separate simulations in which only one of the three factors was varied while the others were kept constant at their 1960s reference.

3 Results

Recent changes in climate have lead to substantial changes in the magnitude and distribution of precipitation and potential evapotranspiration, the two main climate drivers in the Budyko model (Eq. 1) that determine how average precipitation is partitioned between evapotranspiration and streamflow. Average precipitation during the reference period shows a general decrease towards the east (Fig. 2a). Superimposed on this large-scale pattern are local areas with higher precipitation along the coastal areas in the west and/or in mountainous regions. Changes in average precipitation over the study period show a strong north–south gradient (Fig. 2b): most of the Mediterranean, in particular the Iberian Peninsula, shows a decline in precipitation, whereas northern Europe, in particular the British Isles, the Scandinavian peninsula, and Finland, have seen strong increases in average precipitation regionally exceeding 20 %. In contrast to precipitation, potential evapotranspiration shows a strong latitudinal gradient (Fig. 2c) with lower values (PET around 400 mm yr⁻¹) in Scandinavia and higher (regionally exceeding 1000 mm yr⁻¹) in the Mediterranean. Changes in potential evapotranspiration (Fig. 2d) are predominantly positive and highest in central Europe, reflecting the higher increase in average temperatures and shortwave radiation. In general, these strong changes in climate forcing (P and PET) are likely to be reflected in continental-scale patterns of changes in water availability.

Figure 4 Simulated water balance partitioning over the period 1960–2010. Left panels show the mean evapotranspiration and streamflow in 1960 (a and c, respectively), while the right panels indicate the change between the periods 1960 and 2010 for evapotranspiration (b) and streamflow (d).

In addition to climate, land use and land cover in Europe have also seen large-scale shifts over the past 60 years, albeit on a more local scale. Figure 3 shows the mean forest and urban fraction for the reference period, as well as the fractional change over the period 1960–2010. While forest cover is widespread over most of Europe (Fig. 3a), most extensive forest regions can be found in central-western Europe, Sweden, and Finland. Forest cover has increased considerably over most of Europe (Fig. 3b) following abandonment of less-productive agricultural areas and intensification of forestry and forest management, with Sweden (Ericsson et al., 2000) and the Mediterranean region showing the strongest changes. It should be noted that areas where forest cover has declined are virtually absent. This is also true for change in urban areas. The average urban fraction is highest in central-western Europe (Fig. 3c), in particular in Belgium, the German Ruhr area, and the Netherlands. This is also the region that has seen the strongest increase (Fig. 3d). Changes in urban area are generally more localized in nature than changes in forest cover.

Figure 5 Validation of simulated hydrological fluxes across Europe. (a) Simulated ET average over the 10-year periods 1990, 2000, and 2010. (b) Observation-based ET average over the period 1980–2017 from GLEAM version 3.1 (Martens et al., 2017). (c) Simulated changes in streamflow between the periods 1960 and 2000. (d) Observed changes in streamflow over the period 1962–2004 taken from Stahl et al. (2010, their Fig. 2). Note the difference in units between simulations (c) and observations (d) because the approach followed in this study does not allow for normalization by interannual streamflow variability. Observed trends might also be calculated for shorter periods within the period 1962–2004. It should also be noted that ET validation is done for the mean flux, whereas streamflow is validated on the rate of change rather than the mean.

Patterns of mean and changes in evapotranspiration and water yield were calculated by forcing the Budyko model with subsequent 10-year averages of climate forcing and land use at a 1×1 km resolution. Figure 4 shows the resulting continental-scale patterns. The mean evapotranspiration in the reference period (Fig. 4a) is highest in central Europe, locally exceeding 600 mm yr⁻¹, in regions with topographically enhanced precipitation and/or forest cover. The Nordic countries and the Iberian Peninsula generally have lower values (< 400 mm) due to more pronounced energy and water limitation, respectively. Changes in evapotranspiration show a strong latitudinal gradient (Fig. 4b). Changes exceeding +15 % are found in large parts of Scotland, Sweden, Finland, and Estonia, whereas most of central-western Europe shows a smaller increases of the order of 10 %. Decreases of similar magnitude occur in parts of the Iberian Peninsula and Italy. Average streamflow (Fig. 4c) is highest in central-western Europe (locally exceeding 600 mm yr⁻¹), in particular in mountainous areas that receive larger amounts of precipitation. Streamflow of less than 150 mm yr⁻¹ is found in the large parts of Sweden, Finland, Spain, Romania and Bulgaria. Changes in water yield (Fig. 4d) show a roughly similar pattern to changes in evapotranspiration; however, the changes are much stronger in magnitude. Decreases in the Mediterranean locally exceed −45 %, where increases in Sweden and Finland exceed +45 %. Both the changes in evapotranspiration and streamflow show considerable regional variability superimposed on the large-scale patterns.

Figure 6 Two-dimensional quartile distribution of observed versus simulated streamflow. Note that observed streamflow changes are normalized by interannual streamflow variability, whereas simulated changes are normalized by their 1960s values.

In order to assess the quality of the simulated evapotranspiration and streamflow and the changes therein, we evaluate our simulations against observation-based estimates of average evapotranspiration (Martens et al., 2017) over the more recent period 1980–2017 (it should be noted that currently no gridded evapotranspiration estimates are available that cover our complete study period) as well as observed changes in streamflow reported by Stahl et al. (2010) that cover most of our study period. The pattern of simulated ET (Fig. 5a) closely resembles the pattern as produced by GLEAM version 3.2a (Martens et al., 2017, data shown in Fig. 5b). It should be noted that this comparison is added for reference only and should not be seen as a validation: GLEAM is not a strictly observational dataset, and it does not necessarily provide better long-term estimates of ET for forest and urban areas. The Budyko model produces slightly lower values in eastern Europe and the Iberian Peninsula but slightly higher values in Sweden and Finland. At the regional scale, our simulations show more variability due to the higher resolution of the forcing and land use datasets. In addition to matching the pattern of average ET, our approach is also able to reproduce the overall pattern of observed changes in streamflow (Fig. 5c, d). The simulations agree with the observed declines in average streamflow in much of southern and central Europe and increases in the more mountainous, coastal, and/or northern regions. The two-dimensional frequency distribution (Fig. 6) confirms the capability of our approach in reproducing the observed trends in Fig. 5d, with a much higher frequency in the outer quartiles along the diagonal (12 % and 9.7 % compared to 6.25 % expectation) than across the diagonal (4.8 % and 2.8 % of catchments). It should be noted that a higher-order validation on trends is subject to more noise than validation on mean fields, and a perfect match should not be expected, also due to the difference in normalization. Figure 7 shows that our simulations also add information with respect to trends in forcing (P and PET), where PET and to a lesser extent P show a predominant increase over all basins, while observed trends centre around zero change. Overall, the validation shows that our simplified approach is able to capture continental-scale patterns in mean and changes in evapotranspiration and streamflow.

Figure 7 Comparison of median change normalized by the interquartile range (IQR) for observed and simulated streamflow and climate forcing. Normalization was done in order to allow for direct comparison of the changes reported by Stahl et al. (2010), who reported change normalized by interannual streamflow variability and other fluxes with change expressed in percentage.

Figure 8 Distribution of the absolute contribution of climate (P and PET) and land use (LU) changes to changes in evapotranspiration over the period 1960–2010. Colours reflect the relative importance of land use (LU, in magenta or RGB 0, 255, 255), precipitation (P, in cyan or 255, 0, 255), and PET (yellow, 255, 255, 0). Grey indicates no data. Each contribution is inversely scaled between the 2nd and 98th percentiles over Europe to reflect its relative importance. As a result, white (255, 255, 255) indicates locations where all contributions are below their 2nd percentile, and black (0, 0, 0) indicates locations where all contributions are above their 98th percentile. The side panels show the absolute contributions of LU, P, and PET and the net change for the selected regions. (a) Southern Highlands (Scotland), (b) Paris metropolitan area (France), (c) Landes forest region (France), (d) Seville region (Spain), (e) central Sweden, (f) southern Sweden, (g) Styria region (Austria), and (h) Smolyan Province (Bulgaria). Domain averages are listed in Table 2.

Figure 9 Distribution of the absolute contribution of climate (P and PET) and land use (LU) changes to changes in streamflow over the period 1960–2010. See the caption of Fig. 8 for an explanation of the colours. The side panels show the absolute contributions of LU, P, and PET and the net change for the selected regions. (a) Southern Highlands (Scotland), (b) Paris metropolitan area (France), (c) Landes forest region (France), (d) Seville region (Spain), (e) central Sweden, (f) southern Sweden, (g) Styria region (Austria), and (h) Smolyan Province (Bulgaria). Domain averages are listed in Table 2.

In order to understand how changes in fluxes are driven by local changes in climate and land use, Figs. 8 and 9 show how the contribution of the main drivers (precipitation, PET, and land use) to changes in evapotranspiration (Fig. 8) and streamflow (Fig. 9) varies across Europe. This is done by plotting each contribution (as determined from simulations where the other drivers were kept constant) as a separate RGB component, whereby each contribution is rescaled inversely from the 2nd to the 98th percentiles of its spatial distribution over Europe. The resulting colour map thus has a 3-D colour legend. From the distribution of colours, covering most of the possible colours, it can be readily seen that contributions of individual drivers show a strong variability. Magenta indicates that land use-induced changes in evapotranspiration and streamflow are widespread but generally local in character. Yellow colours occur widely in a latitudinal band between 45 and 54^∘ N, indicating that changes in PET have the strongest influence on water balance partitioning in transitional regions, but less so in water-limited and humid northern regions. Finally, the relative impact of precipitation is strongest above 54 and below 45^∘ N. It should be noted that these continental-scale patterns differ from their changes, which are much more uniform (e.g. PET changes in Fig. 2d are fairly homogeneous).

For a more quantitative regional assessment, the subpanels in Figs. 8 and 9 zoom in on several regions. These further illustrate the strong regional divergence in changes in water flux partitioning. In the southern Highlands of Scotland (Figs. 8a, 9a), a strong increase in precipitation has led to a strong net increase in streamflow of +362 mm yr⁻¹, only slightly counteracted by opposing PET and land use (afforestation) effects. Urbanization in the Paris metropolitan area (Figs. 8b, 9b) has reacted to reduced ET (−18 mm yr⁻¹) but combines with increased P into a significant increase in streamflow (+38 mm yr⁻¹). In the Landes forest region (Figs. 8c, 9c), individual effects are small but combine into a strong (−90 mm yr⁻¹) reduction in water yield. ET changes in the Seville region (Figs. 8d, 9d) are moderate (−53 mm yr⁻¹) due to opposing contributions of precipitation decline and afforestation, but these effects combine into a strong reduction on streamflow (−80 mm yr⁻¹). In Sweden, ET changes (Fig. 8e, f) are stronger in the middle of the country, where widespread afforestation and precipitation increase combine (+95 mm yr⁻¹). As a result, increases in streamflow are stronger in the south (+120 mm yr⁻¹), where land use contributions do not reduce the effect of precipitation increase (Fig. 9e, f). In central Austria (Figs. 8g, 9g), PET increases dominate the net ET change (+45 mm yr⁻¹) but combine with precipitation reduction into a strong reduction of water yield (−108 mm yr⁻¹). In the Bulgarian Smolyan Province (Figs. 8h, 9h), contributions combine into a strong ET increase (+80 mm yr⁻¹) but largely cancel out in the net effect on water yield (−19 mm yr⁻¹). The examples highlight the fact that locally, individual changes are often amplified or counteracted by other changes, but because of the water balance constraint this is only true for impacts on either evapotranspiration or streamflow.

When the results are averaged over the continental scale, land use plays a more important role than suggested by Fig. 8. Table 2 lists the Europe-wide changes in evapotranspiration and streamflow as induced by the three main drivers. While changes in ET induced by precipitation are largest when positive and negative contributions are considered separately, the net effect is smaller since decreases in P in southern Europe are largely balanced by increases in the northern parts. As a result, net effects of land use and PET on ET are comparable to those of precipitation (around 40 km³ yr⁻¹ each), with land use having the largest contribution. These contributions correspond to nearly 1300 m³ s⁻¹, the equivalent of the discharge of a large river. The effects on streamflow differ slightly, with P dominating both the positive and net contributions. When zooming in on the near-natural catchments used by Stahl et al. (2010), a different picture is obtained. The contribution of P is less strong, likely because most of the catchments are located in central-western Europe, where precipitation changes have been modest (Fig. 2b) compared to, for instance, Sweden. The net change in ET is mainly driven by land use and PET. For streamflow changes, P is the largest net contributor at around 4 km³ yr⁻¹, but land use contributes significantly with nearly −2 km³ yr⁻¹. For individual large river basins, such as the Rhine basin shown here, the impacts can differ significantly. Rather than precipitation, land use and PET are found to be the main drivers of changes in streamflow over the past decades. The strong sensitivity of streamflow to past land use changes seemingly contradicts the small land use effects under future land use scenarios for this catchment found in previous studies (e.g. Hurkmans et al., 2009).

4 Discussion

Our results on changes in water balance partitioning over Europe are in line with many more local- or regional-scale studies. In some regions, studies have found few to no trends due to dominance of natural variability on change indicators (Hannaford, 2015). For the 6.5 km² Hupsel Brook catchment in the east of the Netherlands, Brauer et al. (2018) reported no significant trend in annual runoff since the mid 1970s. In one of the few studies on long-term in situ observations of ET, Seneviratne et al. (2012) reported no significant trends of annual ET at the Rietholzbach lysimeter in north-eastern Switzerland. These findings are consistent with the results on changes in ET and Q presented in Fig. 4b and d. Other regions have seen negative trends. The decline in water yield in the Ebro River has been attributed to land abandonment (López-Moreno et al., 2011), whereas precipitation decline has been identified as an additional factor in most of the Iberian Peninsula (Lorenzo-Lacruz et al., 2012). In Austria, increased P and PET have been identified as factors driving ET increase (Duethmann and Blöschl, 2018). In Sweden, Jaramillo et al. (2018) found little change in the ratio ET ∕ P in spite of strong increases in P and PET. Also, these findings are consistent with our results. This shows that even using gridded observations contains consistent information for local-scale change analysis.

The modelling approach followed here is simplified in terms of number of model parameters, land use classes, and the parameterization of climate. While the single model parameter w^* correlates with physical land surface properties, it does not have a direct physical meaning (although expressions can be derived linking Budyko parameters to vegetation and climate characteristics; see Gerrits et al., 2009). Therefore w^* might also change with mean climate conditions, the synchronicity between precipitation and potential evapotranspiration, changing snow conditions (Berghuijs et al., 2014), and/or vegetation phenology (Donohue et al., 2007). This could not be investigated due to a lack of observations in southern and northern Europe. It has also been argued that the success of Budyko approaches can be partly explained by the possible adaptation of vegetation to differences in climate seasonality and soil type (Gentine et al., 2012), which would be a strong argument in favour of using such simplified models. We also use a limited number of land use classes. This number is constrained by both the limited availability of accurate estimates of long-term water balance partitioning for different land use types as well as by the limited number of land use classes in the HILDA land use reconstruction. Nonetheless, our simulations capture the most important land use and climate-induced impacts.

(Stahl et al., 2010, and Fig. 5c, d)

Table 2 Climate and land use contributions to changes in evapotranspiration and streamflow over the period 1960–2010. All units in km³ yr⁻¹. For reference, 1 km³ yr⁻¹ corresponds to an average discharge of 32 m³ s⁻¹. The total area with available data is 4 312 807 km².

The lysimeter observations include land use with some of the highest and lowest reported ET rates, making it unlikely that we underestimate the land use-induced variability in ET. Whereas there can be considerable variability in average ET within land use classes, for instance due to vegetation and/or soil type (Haferkorn and Knappe, 2002), this variability is typically small compared to the possible range of ET over all land use classes. The range in w^* values is also consistent, at least qualitatively, with estimates in previous studies (Li et al., 2013; Greve et al., 2014). Our modelling approach did not explicitly consider effects other than atmospheric temperature as climate drivers of ET. For instance, the impacts of rising CO₂ levels on transpiration (Piao et al., 2007) were not considered, although the effects of CO₂ were found to be small compared to effects of forest stand age (Jaramillo et al., 2018). Also, the impacts of agricultural intensification (Liu et al., 2015) and irrigation on ET were not considered, although in some regions the effect of irrigation can be considerable (Siebert and Döll, 2010; Jaramillo and Destouni, 2015). Both can be expected to lead to higher ET and lower Q. While such processes can have strong impacts locally and regionally, other studies have shown small effects under European conditions (e.g. van Roosmalen et al., 2009). It should be mentioned that other, more rigorous, methods have been applied at smaller scales based on multiple working hypotheses (Harrigan et al., 2014) that allow for identification of additional factors driving hydrologic change. The observation that regional disagreement can exist between simulated and observed streamflow changes indicates that more research is needed to fully understand drivers of streamflow change at smaller (regional and catchment) scales.

The model forcing is based on interpolated observations from weather stations. The location of these stations generally follows WMO recommendations (see e.g. Ehinger, 1993), and as a result there is a lack of meteorological observations in, near, or above forests (Frenne and Verheyen, 2016) or in urban areas. Large forest or urban areas, however, are known to impact their own weather, for instance due to enhanced temperature (the well-known urban-heat island effect), cloud formation (as has been observed over the larger French forest regions of Landes and Sologne and cities of Paris and London; see Teuling et al., 2017; Theeuwes et al., 2019), or rainfall (as has been shown by modelling experiments for the Dutch Veluwe forest region; see ter Maat et al., 2013). Such local land cover impacts on climate are unlikely to be represented correctly in the forcing dataset used in this study, which is based on interpolation of weather station data. Also, the quality of the data underlying the E-OBS and HILDA datasets used in this study might differ between countries. As a result, the datasets might induce “jumps” near to borders, as can be seen in some of the maps. These inconsistencies will likely be fixed in future releases of the datasets and do not impact the overall conclusions of this study.

The model forcing of potential evapotranspiration is determined using the Penman–Monteith parameterization (Harris et al., 2014), which accounts for temperature, radiation, humidity, and wind speed effects on evapotranspiration. The benefit of this approach is that it is the most physical model for potential evapotranspiration, but the larger number of variables involved also increases the risk of spurious trends. Routine observations of net radiation, needed to force more complex parameterizations such as the Penman–Monteith equation, are only available for the most recent decades from either stations or satellite. Often, they are calculated from other (uncertain) input data. This raises the question whether decadal trends in radiation (i.e. global dimming and brightening; see Wild, 2016) are correctly represented in long-term PET datasets based on Penman–Monteith. Potentially, trends in PET might be underestimated. A major disadvantage of simpler temperature-based methods is that, while they correctly follow the intra-annual variations in energy, they might be too sensitive to interannual and decadal variations in temperature that are independent of radiation trends (Sheffield et al., 2012). The difference between temperature-based and more physical representations will be minimal, in particular in drier (semi-arid) regions with seasonal water limitation due to the reduced sensitivity of ET to PET (van der Schrier et al., 2011). It has also been reported that temperature-based methods such as Thornthwaite do not always give the strongest increase in PET in a warming climate when compared to other more physically based methods (Prudhomme and Williamson, 2013), suggesting that PET-induced changes in water balance partitioning should be interpreted with care.

Changes in climate and land use generally affect both the average evapotranspiration and streamflow. But whereas changes in evapotranspiration are needed to explain changes in streamflow, the socio-economic impact relates more directly to streamflow since this reflects average freshwater availability. This is of particular relevance in the Mediterranean region, where a decline in water yield or streamflow reflects a decrease in water available for irrigation and agricultural production downstream. Our results indicate that land use changes in the more mountainous areas in the Mediterranean have contributed significantly to reductions in streamflow. Conversely, increasing streamflow in northern Europe might be beneficial to other sectors such as the hydropower industry. The finding that land use change effects are of similar magnitude to climate change effects on water availability also has important implications beyond the yearly average values. Extremes will likely also be impacted by land use, yet current drought projections for Europe (Forzieri et al., 2014; Samaniego et al., 2018) or assessments of changes in floods (e.g. Hall et al., 2014) do not take into account past and/or future land cover changes. Not accounting for land use change will likely lead to regional overestimation or underestimation of changes in water availability. Therefore, land use change impacts on evapotranspiration and streamflow need to be considered in conjunction with climate change impacts.

5 Conclusions

In this study, we investigated the role of changes in land use and climate in Europe from 1960 to 2010 in average evapotranspiration and streamflow. In our modelling approach, we combined a state-of-the-art land use reconstruction with gridded observational datasets of climate forcing and a Budyko model constrained with ET observations from several long-term lysimeter stations. Based on the model results, it was shown that land use changes have had net impacts on evapotranspiration that are generally comparable in size to those caused by changes in precipitation and potential evapotranspiration. Evapotranspiration increased in response to land use (mainly large-scale reforestation and afforestation) and climate change in most of Europe, with the Iberian Peninsula and other small parts of the Mediterranean being exceptions with negative trends. Streamflow changes were dominated by a strong positive contribution of precipitation increases in northern Europe. Land use and potential evapotranspiration had smaller effects of opposite sign, resulting in small net streamflow changes over Europe. The analysis revealed considerable complexity at smaller scales, with most of the possible combinations between positive and negative contributions of precipitation, land use, and potential evapotranspiration occurring at some locations. This was true for effects on evapotranspiration and discharge. Most pressingly, we find that in much of the Mediterranean, land use and climate change combine to further reduce streamflow and water availability.

Appendix A

Figure A1 Location of the stations and sites listed in Table 1. Triangles indicate sites with forest observations. Other symbols are as in Fig. 1.