Blue in green: forestation turns blue water green, mitigating heat at the expense of water availability

All simulations presented in this paper were produced with the Canadian RCM (CRCM5; Šeparović et al 2013, Martynov et al 2013) coupled to the Canadian Land Surface Scheme (CLASS; Verseghy 1991, Verseghy et al 1993). The model setup used is described in greater detail in the sensitivity study (Asselin et al 2022). In this paper, we analyse an ensemble of 12 regional climate simulations (figure 1). These runs are dynamically downscaled global simulations from the Max Planck Institute for Meteorology ESM version 1.2 (MPI-ESM-1.2; Mauritsen 2019). In particular, we used the four first members (rXi1p1f1, X = 1-4) of the lower-resolution ensemble (MPI-ESM-LR-1.2). For each of the four member realisations, we set up three configurations: present GHG and land cover ($G_p L_p$, or control), future GHG but present land cover ($G_f L_p$), and both future GHG and land cover ($G_f L_f$). These three configurations allow us to isolate the respective effects of LUC and GHG on any variable of interest (see appendix A.1 for details).

Zoom In Zoom Out Reset image size

All runs were performed with a horizontal resolution of 0.11^o over the Europe CORDEX domain (https://cordex.org (accessed on 30th April 2024)), except that 30 additional grid points were added to the west side of the computational domain so that all land grid cells are free of spurious spatial spin-up effects (Matte et al 2017). The analysis presented here excludes a 7 year temporal spin-up, such that climatological means are calculated over 1986–2015 (present GHG) or 2071–2100 (future GHG).

This paper focuses on the shared socio-economic pathway SSP1-2.6. In this scenario, global population growth levels off by mid-century, diets shift away from animal products, agricultural yields continue to improve and policies are enacted to halt deforestation, promote afforestation and restore degraded forests, resulting in increased forested area and in substantial reductions in grazing land and agricultural land worldwide (Doelman 2018). In Europe, however, the food crop reduction trend is largely canceled by a massive ramp up of biofuel production (Hurtt 2020). In terms of GHG emissions, SSP1-2.6 includes strong mitigation policies limiting global warming to about 2 °C above pre-industrial levels by the end of the century (Tebaldi 2020).

To quantify the effects of LUC, we need a spatial map of the present land cover as well as the projected changes to that land cover in the future scenario. For the present land cover, we use the version 1.1 of the LANDMATE dataset (Reinhart et al 2022a, 2022b), which provides fractions of 16 plant-functional types (PFTs) over Europe for the year 2015. Hoffmann et al (2023) developed a land-use translator that generates annual land cover maps based on future projections of land use transitions for different SSPs derived from the Land Use Harmonized Dataset Version 2 (LUH2; Hurtt 2020). In this paper, we use their map for the SSP1-2.6 scenario in year 2100 (Hoffmann et al 2022). For consistency with the regional climate simulation ensemble, we chose the 0.1 degree resolution versions of the land cover maps. Note that these maps are publicly available (Hoffmann et al 2022, Reinhart et al 2022a).

To implement these land cover maps in our model, it is necessary to convert the LANDMATE PFTs to our model's PFTs. The conversion is straightforward as the PFT categories of LANDMATE and CRCM5 are very similar (see table SI-LC1). To simplify the analysis, we aggregate the resulting PFTs into six land cover categories: broadleaf trees (btree), needleleaf trees (ntree), shrubs, grasses, crops, and urban. By construction, the land-use translator does not change the fraction of bare soil (Hoffmann et al 2023). This land category was thus discarded from our analysis.

Figure 2 shows how land cover changes between 2015 to 2100 under SSP1-2.6. Consistent with ambitious carbon sequestration efforts, there is a strong forestation signal in Europe. Broadleaf trees account for most of the increase (about 0.2 Mkm²), followed by needleleaf trees (about 0.14 Mkm²). The land cover type showing the most significant decline is grasses (about −0.28 Mkm²), followed by shrubs (−0.09 Mkm²). There are several regions of strong cropland expansion and reduction, but the net area change is rather weak. Finally, urbanization is comparatively negligible on the continent scale ($\lt$0.01 Mkm² ) and will therefore be ignored in the rest of our analysis. In terms of land cover change, the SSP1-2.6 may be thought as a more realistic version of the idealized FOREST experiment of the CORDEX Flagship Pilot Study Land-Use and Climate Across Scales (LUCAS) phase 1 (Davin 2020), where all vegetated area were turned into forests. Most of the positive changes in land cover come from forestation, and most of the land replaced is grassland.

Zoom In Zoom Out Reset image size

In this section, we assess the relative contributions of LUC and GHG emissions to end-of-century summertime heat extremes. To do so, we compute the 95th percentile of the maximum daily temperature (TX95) during the summer months (JJA) across the three distinct model configurations of land cover and GHG concentrations in our simulation ensemble (see detailed methodology in appendix).

Figure 3 shows how TX95 is modified by LUC (left panel), GHG (middle panel), and their combined influence (right panel). LUC cause a widespread alleviation of heat extremes over large parts of central and eastern Europe, where most of changes in land-use occur. The heat mitigation pattern is land-confined and distinctly structured. In fact, some of the fine-scale patterns can be recognized from LUC maps—the thin forestation line stretching across southwestern Russia, for instance—hinting that LUC have highly local effects on heat extremes. Of all LUC categories, $\Delta_\text{LUC}$TX95 is best correlated with broadleaf forestation and shrub reduction (R = −0.57 and R = 0.44, respectively; figure SI-SM1). The cooling effect of LUC is substantial: regions undergoing intense broadleaf forestation respond with a reduction of the TX95 exceeding 1 °C.

Zoom In Zoom Out Reset image size

GHG emissions, on the other hand, cause a widespread aggravation of heat extremes (middle panel). The GHG response is much smoother and extends to all corners of the domain, including over bodies of water. However, the increase in TX95 is mostly concentrated in the semi-arid Mediterranean region and at high latitudes, leaving a region of modest warming (and even localized mild cooling) in Central Europe. As such, when the effects of both LUC and GHG are combined (right panel), one finds that summertime heat extremes are reduced over large swaths of continental Europe. In other words, under SSP1-2.6, LUC (primarily broadleaf forestation) trump GHGs in mitigating summertime heat extremes. Despite GHG-driven global warming, end-of-century heat waves are milder than in the present thanks to forestation.

Using contemporary population data, our estimates suggest that by the end of the century, about half of the European population could experience milder heat waves under this scenario (figure SI-T1; appendix A.3). In contrast, little more than 1% would benefit from such relief in the absence of forestation efforts. Hence forestation could qualitatively alter the experience of climate change for a significant portion of the European population.

What are the mechanisms by which LUC mitigate heat? To tackle this question, we analyse how LUC modify the surface energy balance. Since we are interested in heat extremes, we focus the analysis on hot days, which we define as summer days with a maximum daily temperature (TX) exceeding its 95th percentile (TX95) for a given configuration (see detailed methodology in appendix).

Figure 4 shows how LUC alter the surface energy balance during hot summer days. It is clear at first glance that the magnitude of changes in the latent (LH) and sensible heat (SH) fluxes is far greater than changes in shortwave (SW) and longwave (LW) radiation. In other words, LUC predominantly affect turbulent fluxes, with second-order effects on radiative fluxes. Ground heat fluxes (not shown) are dwarfed by both turbulent and radiative fluxes, and are therefore not considered in our analysis.

Zoom In Zoom Out Reset image size

The spatial patterns of LH increase are approximately matched by decreased SH, implying a strong shift in the partition of turbulent fluxes. The fine-scale spatial structure of this shift bears resemblance to the map of heat extreme mitigation due to LUC (figure 3). To quantify this observation, we calculated correlations between LUC-induced, hot-day energy fluxes and maximum daily temperatures (ΔTX; figure SI-SM2). Patterns of ΔLH match ΔTX very well (R = −0.76), with ΔSH following closely behind (R = 0.65). This means that on hot days, regions where LUC boost evapotranspiration tend to experience stronger heat relief. In effect, LUC redirects solar energy away from SH, which would otherwise result in warming of the atmosphere. This redirected energy is instead utilized to break hydrogen bonds, i.e. evaporate water, without direct temperature change. The net result is the well-known phenomenon of evaporative cooling.

Heat extremes match patterns in the longwave signal even better (R = −0.80) than latent heat fluxes (R = −0.76). This is not surprising as outgoing longwave radiation is slaved to surface temperature. Via Planck's feedback, any change of surface temperature is resisted by the adjustment of upwelling longwave radiation. While the biophysical properties associated with LUC in our model do not directly impact the emission of longwave radiation, LUC reduce surface temperature by boosting evapotranspiration. As a result, upwelling longwave is reduced (net LW increased) and LUC-induced evaporative cooling is damped. That said, the amplitude of ΔLW is weak compared with ΔLH and ΔSH. Despite its strong correlation with heat extremes, longwave radiation only plays the second-order role of damping the first-order effect of evaporative cooling.

Albedo is often the main property of interest when considering the biophysical effects of forestation, since trees tend to be darker than grasses (Hasler 2024). During winter and spring, albedo exerts a strong influence on temperature as snow covers forested and deforested regions (figure SI-T2). However, in the present case of summertime and with LUC dominated by broadleaf forestation of grasslands, albedo is not so important. In our model, broadleaf trees have a marginally smaller visible albedo than grasses (α = 0.05 vs 0.06; see table SI-LC1) , hence LUC induce only weak changes in shortwave radiation. In fact, the correlation between LUC-induced, hot-day ΔSW and ΔTX is negative (R = −0.32), meaning that excess shortwave radiation tends to be co-located with reduced heat extremes. This is the opposite of what would be expected in a situation of shortwave-driven warming.

In summary, the main effect of LUC on the hot-day surface energy balance is to shift the partition of turbulent fluxes away from sensible heat and towards latent heat fluxes. Incident solar energy that would otherwise warm the surface and the air above is instead utilized to evaporate water, hence cooling results. Planck's feedback acts to reduce this cooling effect, but its magnitude is relatively weak. Albebo changes are insufficient to drive appreciable shortwave-driven warming.

We have seen in previous sections that LUC increase evapotranspiration (E) during hot summer days. This remains true for seasonal averages all year round (green bars in figure 5): LUC boost green water fluxes by 0.046 mm d⁻¹ on average over the year. This rather small number conceals important seasonal variability. Since vegetation is most active during the warmer months of the year, there is a strong seasonal cycle in ΔE. While LUC have nearly no effect during winter, they increase E by more than 0.1 mm d⁻¹ on average during the summer. These seasonal averages have weak inter-member variability (see limited spread in the black dots depicting individual member values, figure 5). The highly structured spatial patterns of ΔE (figure 6) and their resemblance to LUC patterns suggest that the evaporative response is dominated by local changes in the biophysical properties of the surface, with secondary dependence on climate variability and local/regional atmospheric feedbacks.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

LUC also increase annual precipitation by 0.032 mm d⁻¹ on average (purple bars in figure 5). The seasonal cycle of ΔP resembles that of ΔE although with weaker amplitude, ranging from negligible wintertime change to peak boost in summer at about 0.07 mm d⁻¹ . Since LUC-induced differences are calculated between couples of simulations with identical atmospheric and oceanic conditions at the boundaries, ΔP can be interpreted as regional precipitation recycling (Eltahir and Bras 1996, Tuinenburg et al 2020). Prevailing winds being westerly, there is little to no agreement in the sign of ΔP upwind (west) of the land mass (figure 6). The most robust component of the signal is the widespread increase of ΔP over continental Europe during the summer. Compared with ΔE, the ΔP signal is much smoother, hinting at the more non-local character of precipitation recycling (Meier et al 2021, te Wierik et al 2021). The spatial variability of continental summertime ΔP is relatively robust across members, although there is more disagreement than for ΔE in some regions. During the other seasons, there is little agreement in the signal (figure SI-WF1). Nevertheless, the spatiotemporal average shows a robust increase in annual precipitation.

Annual LUC-induced ΔE is greater than ΔP by about 0.014 mm d⁻¹ . Annual run-off (R), or the blue water flux, is thus reduced by a comparable amount (blue bars in figure 5). In our model, R is the sum of surface and deep run-off (figure SI-WF2). Surface run-off has a relatively weak contribution to ΔR, except during winter and spring as snow-related dynamics come into effect (figure SI–SN). Total ΔR is mostly driven by changes in deep run-off, which have a less pronounced seasonal cycle. The seasonal cycle of ΔR is out of phase with that of ΔE and ΔP: peak ΔR occurs in spring, with most modest values during summer and fall.

The climatology of soil moisture, M, is significantly influenced by forestation (figure 7). While the northwest portion of Europe experiences only superficial changes, the south and southeast of the continent show both areas of strong replenishment and drying, especially north of the Black Sea and in the Hungarian plains. LUC contribute to the drying of these regions because LUC-induced precipitation recycling does not match the blue-to-green water flux shift. Furthermore, since these regions have relatively dry soils in the first place, the relative changes in soil moisture can be drastic. In particular, about 1% of the significant-LUC grid cells undergo a drying exceeding 25%, and 0.1% dry by 40% or more (figure 7, middle panel).

Zoom In Zoom Out Reset image size

Inspection of the vertical profile of spatially-averaged ΔM, depicted in the right panel of figure 7, is also instructive. While LUC cause no net moisture change in the first top meter of soil, a significant drop in moisture becomes apparent below 1 m. The drop of moisture below 1 m is compelling evidence that forestation modulates the water budget by pumping of water deeper in the soil than grasses and crops would. In our model, broadleaf trees have 2 m deep root whereas grasses, which are the main land cover replaced by forestation, have 1.2 m deep root (table SI-LC1). With broadleaf-tree forestation, soil moisture between 1.2 and 2 m becomes available for pumping. As a result, forestation-induced evapotranspiration drains additional moisture between 1.2 and 2 m. This slows percolation of soil moisture to deeper layers and hence reduced deep drainage (negative ΔR). The LUC-induced moisture profile has an about 2 kg m⁻³ moisture deficit extending from about 2 m down to almost the whole soil column. In summary, moisture profiles reveal a clear signature of root deepening with forestation, and enhanced evapotranspiration by broadleaf trees is associated with reduced deep soil moisture and drainage.

In this paper, we used an ensemble of regional climate simulations to weave together a narrative of the effects of large-scale forestation in Europe under SSP1-2.6. In essence, we found that forestation brings heat relief by turning blue water green, compromising water availability in some already-dry regions. Ambitious forestation plans thus introduce a thorny trade-off between heat mitigation and water availability (Jackson et al 2005). One the one hand, large-scale broadleaf forestation would boost carbon sequestration and green water fluxes, both mitigating the negative impact of extreme summer heat. On the other hand, more green water means less blue water in the soil, streams, reservoirs and aquifers essential to support human life and other ecosystem services. To complicate matters, forestation has non-local effects on precipitation. One region's decision to plant trees for heat relief may affect another region's hydrology.

How well does this narrative hold against other lines of observational and modeling evidence in the literature? Over large parts of continental Europe, we find that LUC trump GHG, resulting in milder heat waves at the end of the century compared with today (figure 3). Using four ESMs, Hirsch et al (2018) also found LUC to contribute significantly to the heat extremes in low-emission scenarios, although typically less so than other forcings combined and with more variability in the sign of change. There is nevertheless considerable observational evidence that the presence (absence) of forests is associated with reduced (increased) heat extremes (Li et al 2015, Alkama and Cescatti 2016). Lejeune et al (2018) used an observationally-constrained climate simulation ensemble to show that historical deforestation increased the intensity of hot days in northern mid-latitudes. Similarly, satellite observations of the recent vegetation greening in China were associated with reduced mean and extreme temperatures (Cao et al 2023). There is further satellite-based (Schwaab et al 2020) and modeling-based (Breil et al 2023) evidence that broadleaf trees provide more heat relief than needleleaf trees, as was found here.

To identify the mechanisms responsible for the heat relief, we analyzed how LUC alter the surface energy balance (section 3.2). Consistent with in-situ and remote sensing observations (Bright et al 2017), we find that LUC primarily affect non-radiative processes (figure 4). Duveiller et al (2018) estimated the effects of different various LUC on the surface energy balance using remote sensing observations. Their deciduous broadleaf forest to grasses transition shows the same sign of change as our analysis for all energy fluxes, except for their weak sensible heat signal. Li et al (2015) also demonstrate that during summer, the effects of enhanced evapotranspiration outweigh those of albedo, leading to lower daily temperatures over northern mid-latitude forests compared to open lands. Likewise, Zeng et al (2017) identify evapotranspiration as the main agent driving biophysical cooling associated with Earth's greening.

In section 4, we explored the implications of large-scale forestation on the water cycle. In our simulation ensemble, forestation boosts both evapotranspiration and precipitation (figures 5 and 6). This aligns with observational studies showing that LUC in Europe over the past decades—mainly forestation of rain-fed agricultural lands—have resulted in increased evapotranspiration (Teuling et al 2019) and precipitation (Meier et al 2021). Similar changes have been observed following afforestation of the Chinese Loess Plateau (Tian et al 2022). Our resulting run-off signal is weaker and patchy, but negative on average over regions of significant LUC, consistent with the well-known drying of streams typically following forestation (Bosch and Hewlett 1982, Jackson et al 2005). We also find that LUC-induced changes in the water cycle can cause significant drying of the soil column in already-dry regions (Ge et al 2023; see figure 7). Hoek van Dijke et al (2022) estimated that global tree restoration could decrease water availability by up to 38% in some regions, which is on par with our findings (figure 7).

While this study provides a self-consistent illustration of the linkage between forestation, climate and hydrology, a myriad of limitations hinders its direct use for policy decisions. First and foremost, our study captures a small fraction of total uncertainty. We explored LUC under SSP1-2.6, which, although a standard scenario in climate modeling (O'Neill 2016), represents just one plausible outcome among countless possible futures. It would be surprising, in fact, if heat extremes were also alleviated in less optimistic scenarios, which include both more emissions and less forestation. Next, using a four-member ensemble of MPI-ESM simulations is a good start for accounting for this ESM's internal variability, but it is definitely insufficient. We chose this ESM in particular because of its plausible 3 °C climate sensitivity, but a single ESM is cannot capture the breadth of plausible outcomes, even within a specific scenario. A few multi-ESM comparison studies have been conducted specifically on the subject of LUC (notably Lawrence 2016, Hirsch 2018, Lejeune et al 2018, Seneviratne 2018), but high intermodel variability makes an in-depth exploration of mechanisms much more challenging.

Furthermore, the global climate data was downscaled using a single RCM. There is also considerable variability in the response to LUC among RCM-LSM (Davin 2020, Asselin et al 2022). Fortunately, the ongoing second phase of CORDEX Flagship Pilot Project LUCAS is expected to provide a more detailed analysis of variability to be expected across RCMs and land surface model using the same land cover maps. Finally, the land surface model used here, CLASS, remains a relatively crude parametric representation of vegetation, which in reality behaves in a highly complex manner. To this we must add the high uncertainty related to the creation of land cover maps, especially their projection in the future (Reinhart et al 2022b, Hoffmann et al 2023).

To estimate the effects of greenhouse gases (GHGs) and land-use changes (LUC) on a variable of interest, we compute how this variable changes between simulation ensembles with different configurations, i.e. land cover and GHG concentration (figure 1). For instance, the effect of GHGs on temperature (T) is given by

$$\begin{align} {\Delta _\text{GHG}} {T} = {T}_{G_f L_p}-{T}_{G_p L_p},\end{align} \tag{ 1 }$$

where the absence of a superscript on G signals that we summed over the N = 4 members to obtain the ensemble average, for instance:

$$\begin{align} {T}_{G_p L_p} = \frac{1}{N}\sum_{X = 1}^N {T}_{G^X_p L_p}.\end{align} \tag{ 2 }$$

Similarly, the temperature change due to LUC reads as

$$\begin{align} {\Delta_\text{LUC}} {T} = {T}_{G_f L_f}-{T}_{G_f L_p}.\end{align} \tag{ 3 }$$

Finally, to obtain the combined influence of LUC and GHGs, we compute:

$$\begin{align} {\Delta_\text{GHG+LUC}} {T} = {T}_{G_f L_f}-{T}_{G_p L_p}.\end{align} \tag{ 4 }$$

This paper focuses on summer heat extremes, which we define using the 95th percentile of daily maximum temperature (TX95). To obtain this quantity, we sort the maximum daily temperatures (TX) of all days of June, July and August (JJA) from all four members, and extract the 95th percentile. This is done independently for each of the three model configurations, i.e. land cover and GHG concentration, such that each configuration is associated with a unique TX95. The effects of GHG and LUC on heat extremes can then be calculated by replacing T with TX95 in equations (1) and (3), the result of which is displayed in figure 3.

To provide a rough estimate of the population exposed to milder heat extremes, we use contemporary population data from the version 3.3 of the History Database of the Global Environment (HYDE; Klein Goldewijk et al 2017) . The heat extreme maps of figure 2 are re-gridded to the 2023 population map using bilinear interpolation (figure SI-T1). For this analysis, only European countries are considered, excluding Russia, Turkey, Georgia, Azerbaijan, and Cyprus. The total 2023 population included in this domain is about 575 million, of which 277 million (i.e. 48%, or about half) would experience milder heat extremes in simulations with both emissions and LUC. When only the effects of emissions are included, this number drops to 8 million or about 1%.

The analysis of section 3.2 focuses on hot summer days. We define hot days as summer (JJA) days with a maximum daily temperature (TX) exceeding its 95th percentile (TX95) for a given configuration. Therefore, once TX95 is known, a list of hot days can be generated:

$$\begin{align} {D}^\text{hot} \equiv \left\{d_1, d_2, \, \ldots d_M \right\}\end{align} \tag{ 5 }$$

where d_i's are all days such that

$$\begin{align} \text{TX}\left(d_i\right) \gt \text{TX95}\end{align} \tag{ 6 }$$

Again, this is done independently for each of the three model configurations. As a result, each configuration is associated with the same number of hot days, M.

Once the hot-day list is generated, the hot-day average of any variable of interest can be obtained. For instance,

$$\begin{align} \text{LH}^\text{hot} = \frac{1}{M}\sum_{t \, \in \, \text{D}^\text{hot}} \text{LH}\left(t\right),\end{align} \tag{ 7 }$$

where LH(t) is the daily average value of latent heat fluxes for a given configuration. By taking the difference between hot-day averages with and without LUC (i.e. following equation (3)), one can estimate the effect of LUC on hot-day latent heat fluxes, $\Delta_{\text{LUC}}$LH^hot (or simply ΔLH in the text of section 3.2), as was done in figure 4.

In this paper, there are several calculations involving horizontal spatial integration (or averaging): the correlations between LUC, climatic and hydrological variables (section 3), the water fluxes budget (section 4.1) and the vertical profile of soil moisture (section 4.2). To emphasize the effects of LUC, all spatial integrals are calculated over areas of significant LUC only. We set the significance threshold to 5%, meaning that a grid point is considered to undergo significant LUC if any of the aggregated land cover categories—broadleaf trees, needleleaf trees, crops, urban, shrubs or grasses—is associated with a change superior to 5% of the tile area.