Vulnerability of European ecosystems to two compound dry and hot summers in 2018 and 2019

In 2018 and 2019, central Europe was affected by two consecutive extreme dry and hot summers (DH18 and DH19). The DH18 event had severe impacts on ecosystems and likely affected vegetation activity in the subsequent year, for example through depletion of carbon reserves or damage from drought. Such legacies from drought and heat stress can further increase vegetation susceptibility to additional hazards. Temporally compound extremes such as DH18 and DH19 can, therefore, result in an amplification of impacts due to preconditioning effects of past disturbance legacies.

Here, we evaluate how these two consecutive extreme summers impacted ecosystems in central Europe and how the vegetation responses to the first compound event (DH18) modulated the impacts of the second (DH19). To quantify changes in vegetation vulnerability to each compound event, we first train a set of statistical models for the period 2001–2017, which are then used to predict the impacts of DH18 and DH19 on enhanced vegetation index (EVI) anomalies from MODIS. These estimates correspond to expected EVI anomalies in DH18 and DH19 based on past sensitivity to climate. Large departures from the predicted values can indicate changes in vulnerability to dry and hot conditions and be used to identify modulating effects by vegetation activity and composition or other environmental factors on observed impacts.

We find two regions in which the impacts of the two compound dry and hot (DH) events were significantly stronger than those expected based on previous climate–vegetation relationships. One region, largely dominated by grasslands and crops, showed much stronger impacts than expected in both DH events due to an amplification of their sensitivity to heat and drought, possibly linked to changing background CO₂ and temperature conditions. A second region, dominated by forests and grasslands, showed browning from DH18 to DH19, even though dry and hot conditions were partly alleviated in 2019. This browning trajectory was mainly explained by the preconditioning role of DH18 on the impacts of DH19 due to interannual legacy effects and possibly by increased susceptibility to biotic disturbances, which are also promoted by warm conditions.

Dry and hot summers are expected to become more frequent in the coming decades, posing a major threat to the stability of European forests. We show that state-of-the-art process-based models could not represent the decline in response to DH19 because they missed the interannual legacy effects from DH18 impacts. These gaps may result in an overestimation of the resilience and stability of temperate ecosystems in future model projections.

1 Introduction

Extreme dry and hot summers in western and central Europe have become more frequent over the past decades (Coumou and Rahmstorf, 2012; Seneviratne et al., 2014), a trend that is expected to continue as global mean temperatures rise (Barriopedro et al., 2011). Hot extremes in Europe are promoted by changes in atmospheric circulation (Coumou et al., 2015; Drouard et al., 2019) and amplified by strong feedbacks between the land surface and the atmosphere, being therefore also associated with severe droughts (Miralles et al., 2014; Samaniego et al., 2018), i.e. compound dry and hot events (DH).

In Europe, DH events usually have strong negative impacts on ecosystems, such as reduced ecosystem productivity (Ciais et al., 2005; Bastos et al., 2020b). After severe drought and heat stress, plant recovery can be lagged, for example due to reduced growth or non-reversible losses in hydraulic conductance or carbon reserve depletion (Ruehr et al., 2019). This in turn may increase vulnerability to another DH if it occurs before complete recovery. Repeated droughts have been linked to increased forest vulnerability in the northern mid-latitudes, albeit with variable responses (Anderegg et al., 2020). Impaired functioning during the recovery period can additionally increase the hazard of subsequent disturbances, e.g. insect outbreaks (Rouault et al., 2006). However, reductions in leaf area, increases in root allocation (McDowell et al., 2008), or reduced growth, caused by reducing evaporative tissue and enhancing water uptake capacity, could also confer an advantage to subsequent droughts (Gessler et al., 2020). It remains unclear whether the increased vulnerability to a subsequent drought can be explained by compounding hazards (e.g. accumulated water deficits or compound heat) or modulating effects due to vegetation responses to the first event.

In Europe, the summer of 2018 was the hottest since 1500 (Sousa et al., 2020) and associated with an unprecedented area affected by drought (Albergel et al., 2019; Bastos et al., 2020a). This DH event resulted in decreases in ecosystem productivity by up to 50 % in central Europe (Bastos et al., 2020a; Buras et al., 2020) and crop yield losses (Beillouin et al., 2020). Part of the central European region affected by the dry and hot summer in 2018 registered another extremely hot and dry summer in 2019 (Boergens et al., 2020; Sousa et al., 2020).

From a hydrometeorological perspective, each of the dry and hot summers in 2018 and 2019 (DH18 and DH19, respectively) can be considered a multivariate compound event in that both high temperatures and strong drought conditions were observed (Zscheischler and Fischer, 2020). Taken together, they can be considered a temporally compound event (Zscheischler et al., 2020). For example, Boergens et al. (2020) have shown that while soil moisture deficits in summer 2019 were not as pronounced as in 2018, total water storage was lower in 2019 due to the water storage deficit resulting from the 2018 event. Given the unprecedented magnitude of DH18, it is likely that at least some ecosystems had not yet fully recovered in 2019. Therefore, from an ecological perspective, DH19 could additionally be considered a preconditioned compound event, where the impact of DH18 may affect the response to DH19 (Fig. 1). Finally, vulnerability to DH events can be further modulated by long-term environmental changes: water savings from reduced stomatal conductance should attenuate drought stress (Peters et al., 2018), but a concurrent decrease in evapotranspirative cooling along with “hotter droughts” may amplify heat stress (Allen et al., 2015; Obermeier et al., 2018).

Figure 1 Conceptual description of the compound DH18 and DH19 events. Dry and hot conditions in both summers were a result of compounding atmospheric drivers (synoptic patterns, preceding climate anomalies, land–atmosphere interactions). The DH18 impacts were preconditioned by seasonal legacy effects in ecosystem functioning from a sunny and warm spring. We hypothesise that legacies from the DH18 event also acted as a precondition of the response to DH19. Further modulating effects (not shown) include impacts of anthropogenic climate change on drivers, hazards, vegetation conditions, and land cover in modulating responses to hazards.

Separating the modulating effects controlled by vegetation responses to global change or by legacies from past disturbances (Kannenberg et al., 2020) and seasonal legacy effects (Buermann et al., 2018) in observations is problematic as it requires considering the compounding effects of multiple drivers (e.g. synergistic effects of heat and drought stress) and separating the role of seasonal and inter-annual legacies both in physical variables (e.g. soil-moisture depletion) and in vegetation vulnerability to those drivers. This can be done by designing counterfactual scenarios to force process-based models, as has recently been done to evaluate seasonal legacy effects of hot and dry springs (Lian et al., 2020; Bastos et al., 2020a). However, it has been argued that Earth system models fail at modelling woody biomass trajectories following droughts (Anderegg et al., 2015), and thus they might miss inter-annual legacy effects from DH events, although no simulations designed to isolate the individual impact of drought over subsequent years have been performed. Alternatively, statistical models can be used to separate such effects based on observational data (Chan et al., 2021).

Using both remote sensing data and an update to the simulations by Bastos et al. (2020a), we attempt to separate these different effects, namely how the combination of hot and dry conditions affected the vulnerability of ecosystems to the two events (multivariate compound event), how the repetition of a dry and hot summer affected the response to DH19 (temporally compound event), and how inter-annual legacy effects due to impacts of DH18 affected ecosystem vulnerability to DH19 (preconditioned compound event). We first use a statistical modelling approach to evaluate whether signs of increased vegetation vulnerability to DH18 and DH19 can be found and to attribute changes in vulnerability to inter-annual legacies and other modulating effects. We then compare observation-based results to updated simulations by state-of-the-art land surface models and dynamic global vegetation models (for simplicity referred to as LSMs) designed to isolate the impacts of DH18 and their legacy effects (Bastos et al., 2020a).

2 Data

2.1 Climate variables

In ecological studies, drought is better characterised by soil moisture anomalies i.e. agricultural drought (Sherriff et al., 2011; Seneviratne et al., 2012; Samaniego et al., 2018), than atmospheric drought indices. We therefore base our drought assessment on two complementary soil moisture datasets. The first is the observation-based soil moisture dataset obtained from SoMo.ml (Sungmin and Orth, 2021), used as reference in this study, and the second, for comparison with SoMo.ml, is given by ERA5 volumetric soil water content (Hersbach et al., 2020).

The SoMo.ml data are generated using a long short-term memory neural network model trained with meteorological forcing from ERA5 and land surface characteristics as inputs and global in situ soil moisture measurements (Dorigo et al., 2011; Zeri et al., 2020) as target variables. The data cover soil moisture in the first 50 cm of the soil and are available at 0.25^∘ lat/long resolution and daily time steps for the period 2000–2019. We remapped the fields to the finer resolution of the MODIS grid and aggregated the data to monthly means. We then subtracted the mean seasonal cycle and long-term linear trend and divided this by the corresponding standard deviation to obtain standardised soil moisture anomalies (SM_anom).

Monthly temperature and volumetric soil water content from the ECMWF ERA5 Reanalysis were obtained from the Copernicus Climate Change Service at 0.25^∘ lat/long resolution (Hersbach et al., 2020) at monthly time steps, selected for the period 2000–2019 (common with SoMo.ml), and remapped to the finer resolution of the MODIS grid using conservative remapping. Standardised anomalies were calculated as described for SM_anom for ERA5 temperature and soil moisture fields (T_anom, ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{ERA}\mathrm{5}}$). Soil moisture anomalies from ERA5 in layers 1–2 (top 28 cm) are used for comparison of drought conditions with those estimated by SoMo.ml, although the two datasets are not fully independent.

2.2 Vegetation and soil data

We used the 16 d enhanced vegetation index (EVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor from the MOD13C1 Climate Modeling Grid (CMG) product. The MOD13C1 CMG provides continuous cloud-free spatial composites from 1 km data projected on a 0.05^∘ lat/long grid (Didan et al., 2015) and was selected for the period 2001–2019. Standardised EVI anomalies (EVI_anom) were calculated following the same approach as for climate variables. The standardisation allows comparing the relative magnitude of anomalies for pixels with distinct temporal variability patterns and with vegetation productivity simulated by LSMs, which have different physical units.

We used land cover distribution in 2018 from the ESA Climate Change Initiative land cover (Kirches et al., 2014) (CCI-LC). The data are originally provided in land cover classes at 300 m spatial resolution and were converted to fractional cover at 0.05^∘ lat/long resolution for forest, grassland, and crop classes using the CCI-LC user tool.

We used isohydricity fields from global satellite measurements from Konings et al. (2017) at 1^∘ lat/long resolution. Anisohydric plants (low isohydricity) show weak regulation of stomatal opening and prioritise carbon assimilation over water savings during droughts. High isohydric plants show strong stomatal regulation of productivity and thereby preserve water at the cost of carbon assimilation during drought.

We use soil available water capacity (AWC) from Ballabio et al. (2016) and Panagos et al. (2012), which used the Land Use and Cover Area frame Statistical survey (LUCAS) topsoil database to map soil properties at continental scale.

2.3 Outputs from land surface and global dynamic vegetation models

Standardised anomalies of gross primary productivity (GPP_anom) and soil moisture (SM_anom) were estimated by the mean of seven land surface models and dynamic global vegetation models (for simplicity referred to as LSMs) between 1979–2019 from an extension of Bastos et al. (2020a) simulations: a baseline simulation for comparison with observations and a factorial simulation to quantify the individual impact of summer 2018 and its legacy effects, when compared to the reference simulation. A detailed description of the models used and the simulation protocol is provided in Appendix A.

First, all model outputs were remapped to a common 0.25^∘ grid, and the multi-model ensemble mean was calculated for the common period with MODIS (2001–2019). The variables were then deseasonalised, detrended, and standardised as was done for the other variables in the study.

3 Methods

3.1 Drought characterisation

We use the observation-based SoMo.ml as a reference dataset to define agricultural drought conditions. Regions with average SM_anom below −1σ (Seneviratne et al., 2012) during summer (June, July, and August, JJA) are considered drought-affected areas during the DH events. Then, a regional domain affected by both DH18 and DH19 events is selected to evaluate the impacts of two consecutive DH events. Within this region most pixels had negative SM_anom and the majority registered ${\mathrm{SM}}_{\mathrm{anom}}<-\mathrm{1.5}\mathit{\sigma }$, but they differ in the magnitude of agricultural drought in DH19. This allows for comparing responses across pixels for different combinations of stress between DH18 and DH19. Since we are interested in evaluating how recovery from DH18 affected impacts of DH19, we limit our analysis to pixels with negative EVI_anom in DH18.

3.2 Compound DH18 and DH19 events

3.2.1 DH18 and DH19 impact characterisation

To characterise different response types to DH18 and DH19, we group pixels using unsupervised clustering of EVI during the two extreme summers. Using an unsupervised method allows for avoiding making assumptions about the magnitude of impacts or the trajectory between DH18 and DH19 (DH18 → DH19) when grouping pixels. For this, we applied a k-means cluster analysis (Hamerly and Elkan, 2003) using two features, corresponding to the EVI_anom fields in DH18 and DH19. Four clusters captured the most significant differences in the impacts of DH18 and corresponding DH18 → DH19 responses: moderate and strong DH18 impacts and moderate and strong impacts by DH19. These clusters were then used to evaluate how LSMs simulate the summer GPP_anom and SM_anom.

3.2.2 Detecting increased vulnerability to drought and heat stress

To better understand the impacts of the two events, we frame them as a combination of temporally and preconditioning compound events (Fig. 1): a sequence of two DH events, whose impacts may be preconditioned by ecosystem vulnerability to DH, especially in the case of DH19. Vulnerability to DH is defined as the impact of the physical hazard (hot and dry conditions) on vegetation and assessed by remotely sensed EVI and modelled GPP anomalies.

The difference between the reference and factorial simulations by LSMs allow for separating the modulating effects of DH18 legacies to the DH19 impacts (dashed arrow in Fig. 1). Separating the legacies in observations is more challenging because the EVI signal at any time step includes signals from both concurrent climate and past legacies and possibly also long-term global change. To do this, we hypothesise that preconditioning effects due to past disturbance legacies (modulating DH19) and global change (modulating DH18 and DH19) should be detectable by changes in ecosystem sensitivity to similar hazards. Increased vulnerability thus corresponds to EVI_anom values lower (more negative or less positive) than those expected for a given drought or temperature anomaly based on past sensitivities. Inversely, increased resistance would result in EVI_anom being less negative or more positive than expected for a given SM_anom.

We assess whether changes in the sensitivity to climate anomalies are detected in DH18 and DH19 using a statistical modelling approach to predict EVI_anom in DH18 and DH19 based on 2001–2017 climate–vegetation relationships. We do this in two steps: first by fitting a linear regression model for mean EVI_anom in each cluster and then, for more detailed insights, by fitting a random forest model at pixel scale, in which we include potential seasonal legacy effects. In both cases, the training period includes other DH events (Ciais et al., 2005; Orth et al., 2016), with similar climate anomalies, particularly 2003, thereby reducing the risk of attempting to predict EVI_anom based on “unseen” climatic conditions.

As a first step, for the spatially averaged variables within each cluster, we fit the following models:

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{\text{EVI}}}_{\mathrm{anom}}^{C_i}=b_{\mathrm{0}}+b_{\mathrm{1}}\times {\stackrel{\mathrm{‾}}{\text{VAR}}}_{\mathrm{anom}}^{C_i},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\text{EVI}}}_{\mathrm{anom}}^{C_i}$ and ${\stackrel{\mathrm{‾}}{\text{VAR}}}_{\mathrm{anom}}^{C_i}$ correspond to the cluster (C_i) spatial average values of EVI_anom and climate variables (growing-season SM_anom or T_anom), respectively. b₀ and b₁ are the coefficients of each linear regression trained on 2001–2017 values. Each model is then used to estimate DH18 and DH19 EVI_anom. Negative model residuals (observations minus predictions) can indicate increased vulnerability, while positive residuals can be a sign of increased resistance.

However, departures from a linear model could also result from non-linear interactions between soil moisture and temperature or from legacy effects from spring (Bastos et al., 2020a; Lian et al., 2020). To account for such effects and evaluate potential spatial asymmetries in the departures from long-term climate–vegetation relationships, we fit a random forest (RF) model using EVI_anom in each pixel (i) from 2001 to 2017 as the target variable, and the corresponding SM_anom and T_anom in spring (March, April, and May, MAM) and in summer (JJA) as predictors:

$$\begin{array}{}\text{(2)}& {\mathrm{EVI}}_{\mathrm{anom}-i}=\text{RF}(T_{\mathrm{anom}-i}^{\mathrm{spr}},{\mathrm{SM}}_{\mathrm{anom}-i}^{\mathrm{spr}},T_{\mathrm{anom}-i}^{\mathrm{sm}},{\mathrm{SM}}_{\mathrm{anom}-i}^{\mathrm{sm}}).\end{array}$$

To reduce the risk of over-fitting due to the small sample size (17 years) and large number of predictors (4), we fit the RF model on 3×3 moving windows centred around each pixel (i.e. 17×9 samples). We assess the model performance outside of the training samples by calculating the out-of-bag scores in addition to the training sample scores. The importance of each predictor is estimated by the Shapley additive explanation values (Lundberg and Lee, 2017). We then predict EVI_anom in DH18 and DH19 using the respective anomalies in $T_{\mathrm{anom}}^{\mathrm{spr}}$, ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{spr}}$, $T_{\mathrm{anom}}^{\mathrm{sm}}$, and ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{sm}}$.

The EVI_anom predicted by the RF model for DH18 and DH19 corresponds to the expected DH impacts from past relationships between the hazards and impacts in Fig. 1. As for the linear case, the difference between the RF model predictions and the actual EVI_anom (model residuals) provides an indication of changes in ecosystem vulnerability to the DH18 and DH19 impacts.

For comparison with LSM simulations, the EVI_anom clusters were remapped to 0.25^∘ by largest area fraction calculation, and subsequently GPP_anom and SM_anom model ensemble means for each cluster were compared with corresponding EVI_anom and ERA5 SM_anom. We first evaluate the linear relationships between the averaged GPP_anom for each cluster and the corresponding climate anomalies for comparison with EVI_anom. Following this, we estimate the legacy effects from DH18 on GPP_anom during 2019 based on the difference between the reference and factorial LSM simulations.

3.2.3 Modulating effects

To understand how land cover can contribute to modulate the impacts of DH18 and DH19, we analyse the land cover composition of each cluster. Given that central Europe is mostly characterised by a very heterogeneous landscape, we calculate land cover selectivity in each cluster for forests, natural grasslands, and croplands. Selectivity is defined as the difference between the probability a given land cover class being present within a cluster compared to its overall presence in the whole region. The probabilities are calculated by fitting a kernel distribution function to the fractional cover fields for the whole region and for separate clusters. Positive (negative) selectivity means that a given land cover type is more (less) likely to be found in a given cluster compared to its overall presence in the region.

Figure 2 Spatial patterns of temperature (T_anom), soil moisture (SM_anom), and EVI (EVI_anom) anomalies during summer 2018 (a–c) and summer 2019 (d–f) for the study region. The study region corresponds to a domain with dry and hot conditions in both 2018 and 2019 summers (DH18 and DH19), delimited by the black rectangle.

For other modulating effects we evaluate how the spatial distribution of EVI_anom residuals for DH18 and DH19 relates to climatic and ecological variables: SM_anom and T_anom in spring and summer, the number of dry months in the year of the DH event and the preceding year (i.e. 2017–2018 for DH18, and 2018–2019 for DH19), EVI_anom in the preceding summer (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$), the number of dry months in a given year and its preceding year (DM), isohydricity (IsoH), and available water capacity (AWC, related to the maximum amount of water available for plants).

We include some of the drivers used to train the temporal climate-driven RF model to diagnose possible changes in the vulnerability explained by stronger vegetation sensitivity to climate anomalies than in the training period. ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ is used to evaluate the preconditioning role of past extreme summers or disturbances (summer is the peak of the growing season in this region). The number of dry months and AWC are also included as they may explain non-linear relationships between SM_anom and vegetation stress. Isohydricity provides a measure of the degree of stomatal regulation by plants. Since many of these variables have strong spatial co-variation (e.g. T_anom and SM_anom, we evaluate their relationships with EVI_anom residuals by calculating the partial rank correlation (Spearman's ρ) between each variable, controlling for the others separately. Since these effects might depend on land cover type, we analyse separately pixels with high and low tree cover.

To further evaluate how inter-annual legacy effects affect long-term vegetation dynamics, we apply a second temporal RF model to pixel-level EVI_anom (Sect. 3.2.2) with ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ as an additional predictor. The model is trained for the period 2002–2017 on 3×3 moving windows and is then used to predict EVI_anom in DH18 and DH19. The resulting model residuals were then compared to those of the climate-driven RF model.

4 Results

4.1 DH18 and DH19 impacts

Following the extreme summer in central Europe in 2018, mild temperatures and strong soil moisture deficits remained until January 2019, when SM_anom returned to normal conditions (Figs. B1 and B2). In central Europe, June 2019 was extremely hot, but July and August 2019 were milder (Fig. B1, Sousa et al., 2020), and soil moisture deficits became very pronounced in July (Fig. B2). In this region, except for during April 2019, the months preceding summer were not particularly dry and were even slightly wetter than average in February, March, and May, the latter also being colder than average. Therefore, the DH18 and DH19 constitute more a sequence of two compound events than a single drought. The areas experiencing severe dry and hot conditions in both summers correspond to a region covering central and eastern Europe and southern Sweden. This region is our study domain and indicated by the rectangle in Fig. 2). Both DH events led to vegetation browning, though negative EVI_anom was more widespread in DH18 than DH19. Within the study region, 79 % of the area showing negative EVI_anom in DH18 (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$) also registered negative EVI_anom in DH19 (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$).

Figure 3 Classification of impact groups within the study region in central Europe. Panel (a) shows the spatial distribution of the four clusters from unsupervised classification of (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) values. The corresponding (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) distribution in each cluster is indicated in panel (b) (circles indicate the spatial mean and lines the spatial standard deviation within each cluster). The corresponding distributions of SM_anom and T_anom pairs are shown in panels (c) and (d), respectively. The grey line indicates similar anomalies in the two DH events. Only pixels with negative ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$ are considered.

The spatial distribution of the clusters resulting from the unsupervised classification based on (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) pairs and corresponding centroids are shown in Fig. 3a and b, as well as the corresponding (${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) and ($T_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, $T_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) (Fig. 3c and d). The four clusters aggregate pixels according to different impacts in DH18 and DH19. One cluster, covering 20 % of the area, includes pixels with moderate impacts in DH18 and further browning in DH19, being therefore referred to as (C_Decline) (dark brown, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$ below the 1:1 line in Fig. 3a). This cluster is associated with mixed cover of forests (10 %–40 %, dominated by needle-leaved forest) and grasslands (15 %–60 %), (Fig. B3). Cluster C_HighV (high vulnerability, covering 15 % of the area) corresponds to pixels experiencing strong impacts in both events and is associated with high grassland and cropland fractions and low tree cover. Pixels with strong impacts in DH18 and weakly negative ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$, i.e. partial recovery in DH19 (C_PRecov, 21 % of the area), are mainly dominated by croplands. Finally, a group of pixels shows moderate ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$ and positive ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$ (C_Greening, 44 %), corresponding mostly to mixed forest–grassland pixels (30 %–65 % of forest, dominated by needle-leaved forest).

All clusters align along proportional DH18 : DH19 values of SM_anom and T_anom, with predominantly negative SM_anom and positive T_anom in both DH events but alleviation of soil moisture deficits and heat stress in DH19 compared to DH18 (Fig. 3). The two recovery clusters (C_PRecov and C_Greening) correspond to pixels with less severe drought conditions and milder temperatures in DH19, and C_Greening corresponds to pixels where dry and hot conditions in DH18 were also more moderate. C_HighV corresponds to pixels experiencing drier and hotter anomalies in both summers and accordingly shows stronger impacts. Cluster C_Decline, however, shows increasing browning in DH19 in spite of drought and heat stress alleviation (Fig. 3). The distributions of climate anomalies for each cluster overlap each other and, in some cases, the 1:1 line, indicating that the intensity of the hazards (temperature, drought) cannot account for the resulting impacts alone.

4.2 Ecosystem vulnerability to DH18 and DH19

All clusters show significant positive linear relationships between summer EVI_anom and SM_anom and negative linear relationships with T_anom in 2001–2017 (Fig. 4). The relationships include the two extreme summers of 2003 and 2015, which had comparable T_anom and SM_anom to DH18 and DH19 in most clusters. However, the long-term sensitivities estimated are robust even if these summers are excluded.

The results correspond to a general summer water-limited regime, especially in clusters C_Decline, C_HighV, and C_PRecov, which show stronger sensitivities to T_anom and SM_anom (slopes in Fig. 4), and higher variance explained by both models (R² 0.58–0.68 for SM_anom and 0.49–0.55 for T_anom). For these clusters, EVI_anom is below the 95 % confidence interval of the long-term linear relationships for DH18 (C_PRecov and C_HighV) and DH19 (C_Decline and C_HighV). SM_anom and T_anom in DH18 and DH19 are generally similar to those of 2003, but DH18 was drier than 2003 in C_PRecov and C_HighV.

Figure 4 Departure of EVI_anom in DH18 and DH19 from long-term climate-driven variability. Relationship between EVI_anom and SM_anom (top row) and between EVI_anom and T_anom (bottom row) for each individual summer between 2001 and 2019 over the study region. The results are shown separately for the four clusters defined in Fig. 3. The black line and shaded areas show the relationship and respective 95 % confidence intervals obtained by ordinary least-squares linear regression between EVI_anom and the respective climate variable for all years between 2001–2017. Values of (EVI_anom, SM_anom) that deviate from the long-term relationships show increased sensitivity to climate anomalies, which can be a sign of increased vulnerability or decline. The colours indicate individual years, ranging from 2001 (red) to 2019 (purple), and square markers indicate 2018 and 2019.

These departures may be related with seasonal legacy effects from the warm spring in DH18 and/or the onset of non-linear responses to heat and drought. To account for these modulating effects, we model long-term (2001–2017) EVI_anom–climate relationships using spring and summer SM_anom and T_anom as predictors using random forest regression (see Sect. 3.2.2). The model is able to predict 48 %–90 % (median and maximum out of bag score across pixels) of the pixel-level temporal variability of summer EVI_anom in 2001–2017 (Fig. B4). Analysis of the variable importance shows that the model estimates summer water limitation and negative legacy effects from spring warming (Fig. B5), consistent with a summer water-limited regime and process-based modelling studies (Bastos et al., 2020a; Lian et al., 2020).

As in the linear case, the RF model estimates less negative or more positive EVI_anom in DH18 and DH19 than observations (Fig. 5). The residuals are below the range of the training period for the high-impact clusters: C_Decline and C_PRecov in DH19 and DH18, respectively, and C_HighV in both (Fig. 5c). In C_Greening, residuals are predominantly positive (i.e. observed EVI_anom more positive than predicted), but still partly overlap with the range of residuals in the training period (Fig. 5).

Figure 5 Spatial distribution of EVI_anom residuals in DH18 (a) and DH19 (b) estimated by the temporal RF model trained for 2001–2017 with spring and summer SM_anom and T_anom as predictors. The corresponding distribution per cluster for each DH event is shown by the boxplots in panel (c). The shaded grey envelope indicates the range of residuals in the training period.

Figure 6 Spatial partial correlation (Spearman) between EVI_anom residuals and environmental variables in DH18 (a, b) and DH19 (c, d) for pixels with high (dark green, top 5 % tree cover fraction) and low (light green, lowest 5 % tree cover fraction) tree cover (a, c). High tree cover (TC) pixels have tree cover fractions above 58 % and low TC pixels have virtually no trees (TC<0.4 %). The variables considered are spring and summer T_anom and SM_anom (indicated by superscripts “spr” and “sm”, respectively), EVI_anom in the previous growing season (EVI_yr−1), plant isohydricity (IsoH), and the number of dry months (DM). Because of the large number of pixels considered, all correlations are significant ($p-\mathrm{val}\ll \mathrm{0.01}$). Panels (b) and (d) show the distribution of residuals for pixels with high and low tree cover.

Pixels with high tree cover tend to show less negative or more positive residuals than pixels with low tree cover (<0.4 %) in both DH events (Fig. 6), but in DH19 the range of residuals in high tree cover pixels is larger than DH18, including pixels with strongly negative values. The partial rank correlation of the spatial distribution of EVI_anom residuals with respect to different explanatory variables is shown for pixels with high and low tree cover in Fig. 6. Given the large number of pixels, all correlations are significant.

In DH18, T_anom in spring ($T_{\mathrm{anom}}^{\mathrm{spr}}$, + for high and low tree cover) and summer SM_anom (${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{sm}}$, − for high tree cover and + for low tree cover) show the strongest relationships with EVI_anom residuals. In DH19, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ (+), $T_{\mathrm{anom}}^{\mathrm{spr}}$, and $T_{\mathrm{anom}}^{\mathrm{sm}}$ (−) show strong correlations, with consistent sign for both high and low tree cover pixels. DH19 residuals of pixels with high tree cover show strong correlation with SM_anom with opposite signs in spring (+) and summer (−) and with AWC (−). In DH19, pixels with low tree cover show negative correlation between IsoH and EVI_anom residuals.

To test whether the importance of ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ is particular to DH19 or if it reflects long-term inter-annual legacy effects of anomalies in vegetation activity, we fit a second temporal RF model where ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ is used as an additional predictor (Figs. B4 and B6). Including vegetation condition in the previous summer improves the predictive power of the long-term RF model (72 %–97 % out-of-bag score, compared to 48 %–90 % for the model trained with climate drivers only). Even though the residuals for the training period are considerably reduced relative to the climate-driven model, the residuals for DH18 and DH19 are comparable.

4.3 DH18 and DH19 impacts simulated by LSMs

The GPP from the LSM multi-model ensemble mean matches the differences in impacts between clusters in DH18 well (Fig. 7a and b) and the temporal evolution of GPP anomalies during the 2018 growing season (April to September, Table 1), with correlations with EVI_anom of 0.74–0.90. Even though the root-mean-squared error (RMSE) is comparable in the two growing seasons, the correlations of GPP_anom with growing-season EVI_anom are much lower in DH19 (−0.09 to 0.43). GPP_anom by LSMs is above average in spring and early summer 2019 for all clusters, and anomalies in DH19 are either more positive or less negative compared to EVI_anom.

LSMs simulate a stronger attenuation of drought compared to the observation-based SM_anom, albeit with consistent relative differences in SM_anom between clusters (compare Figs. B7 and 3). LSMs simulate the temporal evolution of SM_anom well in the two growing seasons, with high correlation with both SoMo.ml and ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{ERA}\mathrm{5}}$ (correlations of 0.81–0.98). The RMSE for simulated SM_anom is generally lower than that of GPP_anom.

Figure 7 Observed and process-based model simulations of 2018/19 impacts. Seasonal evolution of EVI_anom (a) and standardised GPP anomalies (GPP_anom, b) over the 2-year period for each cluster (defined in Fig. 3 and shown for LSM grid in Fig. B7). Panel (c) shows the difference between the reference and factorial simulations and indicates the impacts of DH18 on GPP_anom simulated by models during the event and in the subsequent months until December 2019.

Table 1 Correlation between growing season (gs, April–September) SM_anom simulated by LSMs with SM_anom from SoMo.ml and ERA5 and of EVI_anom with GPP simulated by LSMs.

The sensitivity of GPP_anom to simulated SM_anom and to T_anom (Fig. B8) is consistent with that of EVI_anom in all clusters (Fig. 4), although for C_PRecov and C_Greening the LSMs estimate non-significant negative relationships between GPP_anom and T_anom. The deviations of GPP_anom from the linear response for C_HighV and C_PRecov in DH18 are correctly captured by LSMs but this is not the case for DH19 in C_Decline.

5 Discussion

5.1 Early signs of increased vulnerability

For three clusters covering 56 % of the pixels negatively impacted by DH18, the extremely low EVI_anom in response to DH18 and DH19 could not be predicted from EVI–climate relationships in 2001–2017 (Figs. 4 and 5). These departures reveal increased sensitivity to dry and hot conditions and can be a sign of increased ecosystem vulnerability to such events. However, it should be noted that while we focused on pixels that were negatively impacted by DH18, some pixels in the regional domain selected showed greening, even in DH18 (Fig. 2). These regional asymmetries result in partial regional compensation of the DH18 impacts, as shown in Bastos et al. (2020b).

In both DH18 and DH19, higher tree cover fraction is associated with more positive or less negative residuals (Fig. 6), indicating that trees were more resistant to DH than grasses and crops. The predominance of crops and grasslands in C_HighV, which had strong negative residuals in both events, and of high tree cover in C_Greening also support this effect. Trees can better cope with drought with their deeper rooting depth (Fan et al., 2017) and through the use of carbon reserves to support activity under stress conditions (Wiley, 2020). Moreover, some trees and grasses with stronger stomatal regulation can buffer the drought progression and its impacts by avoiding hydraulic failure (McDowell et al., 2020; Teuling et al., 2010). This is reflected in the small but positive relationship between isohydricity and EVI_anom residuals in pixels with high tree cover.

Increased vulnerability may be explained by modulating effects of global change on vegetation condition (e.g. “hotter droughts”, Allen et al., 2015, Fig. 1) and, in the case of DH19, it may be further linked to inter-annual legacies from the impact of DH18. The first should be expressed by relationships between EVI_anom residuals and climatic variables. The latter are more difficult to assess without comprehensive data about different competing factors, e.g. defoliation or damage from embolism (Ruehr et al., 2019), higher susceptibility to diseases and pests due to reduced health (McDowell et al., 2020) or increased hazard of insect disturbances due to warm conditions (Rouault et al., 2006; Suzuki et al., 2014). The relationships between EVI_anom residuals and ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ provide an approximation but do not allow the identification of the underlying drivers.

In DH18, we find a positive effect of spring warming in vegetation growth, leading to weaker departures from long-term vegetation–climate relationships (observed EVI_anom more positive or less negative than modelled), but with associated water depletion amplifying the impacts of DH18 in summer in pixels dominated by grasslands and crops (low tree cover in Fig. 6). These results are in line with Bastos et al. (2020a) that showed contrasting seasonal legacy effects of warm springs in cropland- versus forest-dominated regions.

On the contrary, spring and summer $T_{\mathrm{anom}}^{\mathrm{sm}}$ in 2019 (or cooling, see Fig. B1) are negative correlated with EVI_anom residuals in both high and low tree cover pixels. This indicates increasing damage from heat stress, for example due to reductions in evapotranspirative cooling (Obermeier et al., 2018) or cascading impacts of compound heat and drought, such as insect attacks (Rouault et al., 2006).

Including ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ in the long-term RF regression model improves the predictive skill for 2001–2017 but does not reduce the residuals in DH18 and DH19. The high correlation between EVI_anom residuals and ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ in DH19 can indicate either that pixels strongly impacted by DH18 were associated with amplified impacts by DH19 (negative residuals) or that pixels affected moderately by DH18 (less negative ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$) were associated with positive residuals, i.e. stronger recovery. Damage to roots and tissues or depletion of carbon reserves from DH18 leading to higher vulnerability to DH19 could explain the positive correlation in high tree cover pixels in C_Decline. Conversely, the moderate DH18 impacts in C_Greening may have resulted in increased resistance to DH19. The strong correlation found in low tree cover pixels is surprising though, as European crop species tend to be annual plants, and annual species can also be found in many grasslands. For these pixels, it is more likely that the positive correlation is explained by management practices, e.g. through earlier harvest or active reduction of stand density in DH19 (Bodner et al., 2015).

C_Decline stands out from the other clusters, in that browning is found in spite of drought alleviation in DH19. The strong negative correlation of residuals with ${\mathrm{SM}}_{\mathrm{anom}}^{\mathrm{sm}}$ and AWC in forest-dominated pixels is counter-intuitive and suggests that other environmental effects not considered in our analysis may modulate DH19 impacts. Insect outbreaks are a potential candidate to explain such effects: the stronger correlation of residuals with ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ in DH19 could reflect increased susceptibility of impaired trees, combined with favourable climatic conditions for insect growth, reflected in stronger negative effects of $T_{\mathrm{anom}}^{\mathrm{sm}}$ in DH19 in high tree cover pixels.

Results from field inventories and forest plots support this hypothesis. Increased tree mortality and insect outbreaks in central Europe during 2018 have been reported (Schuldt et al., 2020). A recent assessment by the German Federal Ministry for Food and Agriculture (BMEL, 2020) reported crown damage in 36 % of all tree types in summer 2019, a 7 % increase compared to 2018 and predominating in trees over 60 years of age. According to this report, the mortality rate in both needle-leaved and broad-leaved trees almost tripled from 2018 to 2019. Although no large-scale data on insect outbreaks are currently available, local authorities in regions where C_Decline is prevalent report an increase in tree mortality from bark beetle infestations: the Environment Ministry of North Rhine-Westphalia in western Germany reported soaring rates of spruce affected by severe bark beetle infestations, from about 1 % in 2018 to over 12 % in 2019 (MULNV-NRW, 2019). In the Czech Republic, rates of spruce damaged by bark beetles more than tripled, leading to increased mortality (Hlásny et al., 2021). In Belgium, a “bark beetle task force” was created in September 2018 by the economic office of Wallonia (OEW, 2018). Increased tree mortality and bark beetle infestations have also been reported in eastern France (ONF, 2020).

5.2 Implications for Earth system modelling

Temperate ecosystems are an important global sink of CO₂ (Pan et al., 2011) and are not usually considered hotspots of drought risk and environmental degradation under climate change (Vicente-Serrano et al., 2020). Our results show that the past two extreme summers in central Europe reveal the first signs of large-scale enhanced vulnerability in response to DH events (C_HighV, C_PRecov) and of potential degradation trajectories induced by consecutive events (C_Decline). Even though it is limited to 20 % of the study area, the patterns in C_Decline highlight the risks associated with more frequent and intense droughts and heatwaves expected in the coming decades (Barriopedro et al., 2011; Boergens et al., 2020; Hari et al., 2020). At the same time, progressive warming conditions can increase the likelihood of compound occurrence of multiple disturbances, such as droughts and insect outbreaks, which are both promoted by warm and dry conditions. Interactions between compounding disturbances can further contribute to forest C losses (Seidl et al., 2017; Kleinman et al., 2019). To anticipate such impacts, process-based modelling of ecosystem response to such events is needed.

The LSMs perform well in simulating the magnitude and evolution of productivity anomalies in 2018 but not in 2019. The recovery simulated by LSMs in DH19 can be partly explained by a strong recovery of modelled soil moisture (Fig. B7) but may also result from limited ability of LSMs to simulate changes in ecosystem vulnerability during the two DH events. The latter is supported by the fact that simulated SM_anom shows good agreement in the temporal evolution of soil moisture anomalies with both observation-based datasets but not of GPP_anom (Table 1).

The comparison of the reference and factorial simulations allows showing that the poor performance in 2019 may be related with interannual legacy effects. LSMs estimate legacies from DH18 only in the early growing season (March to May 2019) but do not estimate any legacy effects in summer (Fig. 7 bottom panel). The poor relationships between EVI_anom and simulated GPP_anom in response to DH19 indicate that processes controlling legacy effects such as damage from embolism, carbon starvation, and resulting tree mortality or disturbances induced by drought and heat, such as insect outbreaks, currently missing in LSMs, likely explain the amplified impacts of DH19.

LSMs are known to have limited ability to simulate drought-induced stress and tree mortality (Wang et al., 2012) and lack impacts of biotic disturbances, although rudimentary approaches have been attempted (Kautz et al., 2018). These model shortcomings add to limitations in simulating soil moisture variability and transitions between energy-limited and water-limited regimes.

6 Conclusions

The summers of 2018 and 2019 were both exceptionally hot and dry over central Europe, and both were associated with widespread vegetation browning and tree mortality events. Here we propose an approach that analyses this event as a combination of three types of compound events (Zscheischler et al., 2020) that consider (i) the compound effects of hot and dry conditions, (ii) the effect of repeated stress conditions in 2019, and (iii) the legacy effects from DH18 impacts in preconditioning the impacts of DH19. Using statistical and process-based modelling, we quantify these effects and identify modulating effects, e.g. land cover composition. This approach can be extended to other types of events that may not fall in a single type of compound event.

Based on remote sensing data, we find signs of degradation trajectories in 20 % of the study area, with vegetation browning in spite of drought alleviation in DH19. We showed that inter-annual legacies from DH18 played an important preconditioning role in amplifying the impacts of DH19. While LSMs simulated the impacts of the first event (DH18) well, they showed limited skill in simulating the impacts of the subsequent compound event (DH19).

Our results show that compounding effects of multiple and repeated stressors and ecological dynamics can result in non-linear and unexpected impacts (Schuldt et al., 2020) that LSMs still cannot realistically simulate. Attribution of inter-annual legacy effects from DH18 and of LSM errors to internal processes (e.g. drought-induced damage and mortality) or others such as insect outbreaks remains challenging because up-to-date datasets on tree mortality and tree carbon reserves or spatially explicit information on biotic disturbances are very limited.

Since extreme DH events are projected to become more common in the coming decades, better understanding the interactions and feedbacks between climate extremes, natural disturbances, and ecosystem dynamics is fundamental to anticipate threats to the stability of forests in the temperate regions and elsewhere. Overlooking these effects may result in an overestimation of the resilience of the CO₂ sink to climate change.

Appendix A: Supplementary methods

Land surface and global dynamic vegetation model simulations

We have used output of gross primary productivity (GPP) and simulated soil moisture from seven models that followed the protocol and extended the simulations in Bastos et al. (2020a) up to 2019. These models are ISBA-CTRIP (Joetzjer et al., 2015), JSBACH (Mauritsen et al., 2018), LPJ-GUESS (Smith et al., 2014), LPX-Bern (Lienert and Joos, 2018), OCN (Zaehle et al., 2010), ORCHIDEE (Krinner et al., 2005), and SDGVM (Walker et al., 2017).

The model simulations were run for most models at 0.25^∘ spatial resolution for the European domain (32–75^∘ N and −11–65^∘ E), following a spin-up to equilibrate carbon pools. For the reference simulation, the models were forced with observed CO₂ concentration from NOAA/ESRL and changing climate between 1979 and 2019 from ERA5 and fixed land cover map from 2010 from LUH2v2 (Hurtt et al., 2011). An additional simulation was run where the models were forced with changing climate, except June–August 2018, where climatological summer climate conditions were used to force the models as described in Bastos et al. (2020a). This simulation, extended up to December 2019, allows for evaluating the direct impact of DH18 and its inter-annual legacy effects.

For more details about the simulation protocol, we refer the reader to Bastos et al. (2020a).

Appendix B: Supplementary Figures

Figure B1 Monthly temperature anomalies during 2018 and 2019. The rectangle indicates the study region.

Figure B2 Monthly soil moisture anomalies during 2018 and 2019. The rectangle indicates the study region, i.e. the areas experiencing drought conditions (${\mathrm{SM}}_{\mathrm{anom}}<-\mathrm{1}\mathit{\sigma }$) during both DH18 and DH19.

Figure B3 Selectivity of different land cover composition for each cluster (Fig. 3). Selectivity is evaluated as the difference between the probability distribution of a given land cover type (forest, a; grassland, b; cropland, c) and the probability distribution of that land cover type in the selected region. If selectivity is positive, the cluster is preferentially composed by the given land cover type and the opposite for negative values. The 2018 land cover classification maps from ESA CCI-LC are used.

Figure B4 Performance of the temporal RF model in predicting EVI_anom given by the out-of-bag scores. Panels (a) and (b) show the scores for the climate-driven RF model, and panels (c) and (d) show the corresponding results for the same model but including ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$ as an additional predictor.

Figure B5 Importance of the four predictors used in the RF model to predict EVI_anom in spring (a, c) and summer (b, d), SM_anom (a, b), and T_anom (c, d) calculated from the Shapley additive explanation values.

Figure B6 The same as in Fig. 5c but for the RF model trained using spring and summer SM_anom and T_anom as predictors, as well as ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{yr}-\mathrm{1}}$.

Figure B7 Panel (a) shows the spatial distribution of the four clusters from unsupervised classification of (${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{EVI}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) values remapped to the coarser grid of LSMs. The corresponding (${\mathrm{GPP}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{18}}$, ${\mathrm{GPP}}_{\mathrm{anom}}^{\mathrm{DH}\mathrm{19}}$) values simulated by the multi-model mean in each cluster are indicated in (b) (circles indicate the spatial mean and the lines spatial standard deviation within each cluster). The corresponding distribution of simulated SM_anom pairs in each cluster are shown in (c). The grey line indicates similar anomalies in the two DH events.

Figure B8 The same as Fig. 4 but for GPP and soil moisture anomalies simulated by a subset of land surface models from (Bastos et al., 2020a) extended up to December 2019.