Identification of compound drought and heatwave events on a daily scale and across four seasons

Compound drought and heatwave (CDHW) events can result in intensified damage to ecosystems, economies, and societies, especially on a warming planet. Although it has been reported that CDHW events in the winter season can also affect insects, birds, and the occurrence of wildfires, the literature generally focuses exclusively on the summer season. Moreover, the coarse temporal resolution of droughts as determined on a monthly scale may hamper the precise identification of the start and/or end dates of CDHW events. Therefore, we propose a method to identify CDHW events on a daily scale that is applicable across the four seasons. More specifically, we use standardized indices calculated on a daily scale to identify four types of compound events in a systematic way. Based on the hypothesis that droughts or heatwaves should be statistically extreme and independent, we remove minor dry or warm spells and merge mutually dependent ones. To demonstrate our method, we make use of 120 years of daily precipitation and temperature information observed at Uccle, Brussels-Capital Region, Belgium. Our method yields more precise start and end dates for droughts and heatwaves than those that can be obtained with a classical approach acting on a monthly scale, thereby allowing for a better identification of CDHW events. Consistent with existing literature, we find an increase in the number of days in CDHW events at Uccle, mainly due to the increasing frequency of heatwaves. Our results also reveal a seasonality in CDHW events, as droughts and heatwaves are negatively dependent on one another in the winter season at Uccle, whereas they are positively dependent on one another in the other seasons. Overall, the method proposed in this study is shown to be robust and displays potential for exploring how year-round CDHW events influence ecosystems.

1 Introduction

Compound events are defined as a combination of multiple drivers and/or hazards that could lead to extreme ecological and socioeconomic impacts of a larger magnitude than the sum of those of the individual components (Zscheischler et al., 2018; Masson-Delmotte et al., 2021). Under the challenge of climate change, extremes, especially those related to temperature, are being observed at an increasing pace (Perkins et al., 2012; Sharma and Mujumdar, 2017; Mukherjee and Mishra, 2021; Kruger and Sekele, 2013) and are expected to become even more severe in the future (Byrne, 2021; Wu et al., 2021). It is critical to study heatwaves concurrently with droughts due to the intensification of negative impacts, such as exacerbating water shortage (Coffel et al., 2019), crop failure and gross primary production reduction (Ciais et al., 2005; Mishra et al., 2020; Zampieri et al., 2017), and wildfire and tree mortality (Libonati et al., 2022; Brando et al., 2014; Reichstein et al., 2013).

Recently, the increasing frequency of compound drought and heatwave (CDHW) events (or compound hot and dry events) has received attention, both globally (Mukherjee and Mishra, 2021; Zhang et al., 2022) and regionally, e.g., in the USA (Mazdiyasni and AghaKouchak, 2015; Alizadeh et al., 2020), India (Sharma and Mujumdar, 2017), Europe (Manning et al., 2019), China (Kong et al., 2020), and Southeast Brazil (Geirinhas et al., 2021).

The identification of CDHW events lies at the core of temporal trend and frequency analyses. However, the coarse temporal resolution and scale inconsistency of the indices used hamper the reliability of existing studies. Droughts are commonly identified on the basis of monthly magnitude variations in historical climate variables (McKee et al., 1993; Ridder et al., 2020; Salvador et al., 2020), and therefore completely ignore the intra-monthly distribution, whereas heatwaves are usually obtained on a daily scale. For example, the 3-month standardized precipitation index (SPI) for meteorological droughts and the daily temperature index for heatwaves are a popular pair of indices used to identify compound events (see, e.g., Geirinhas et al., 2021). This coarse scale of the drought index and the mismatch of scales between the drought and heatwave indices could result in some bias in the actual start and/or end dates as well as the severity of droughts, further affecting the identification of CDHW events. Recent studies have upscaled the timescale of drought events to a weekly, 5 d, or daily scale (Mukherjee and Mishra, 2021; Wang et al., 2020; Li et al., 2020; Mo and Lettenmaier, 2015; Yuan et al., 2023). For example, Mukherjee and Mishra (2021) identified drought weeks as weeks in which the self-calibrated Palmer drought severity index (Wells et al., 2004) falls below the 10th percentile, thereby regarding short dry spells (a week) as drought events. Moreover, Wang et al. (2020) calculated daily standardized precipitation evapotranspiration index (SPEI) values and identified a drought event as a period of consecutive days with an SPEI below a given threshold (such as −1.0). While these identification methods enhance the temporal resolution, they also introduce a new challenge. They may result in a large number of short dry spells and/or mutually dependent ones. Consequently, the droughts identified might not be extreme or independent. Therefore, additional efforts are needed to address the issues of coarse resolution and scale inconsistency.

Most studies on CDHW events exclusively focus on the summer or extended summer season, mainly because heatwaves that lead to deadly heat stress for human beings and ecosystems are only expected in this period of the year (McKechnie and Wolf, 2010; Ummenhofer and Meehl, 2017; Stillman, 2019; Soroye et al., 2020). Although overlooked in the study of CDHW events, ecologists have studied the (ecological) impact of relatively high temperatures (and droughts) in non-summer seasons. For example, warmer winters and earlier springtime drying of soils and forest fuels (fallen leaves) are probably linked to increasing numbers of large wildfires and total area burned in forests in the western USA (Westerling et al., 2006; Littell et al., 2009; Williams et al., 2010), which could weaken the resistance of surviving trees to bark beetle attack (Raffa et al., 2008). Laboratory experiments conducted by Radchuk et al. (2013) and Oliver et al. (2015) and a long-term countrywide data set analysis by McDermott Long et al. (2017) all found that elevated temperatures during the overwintering period could increase the mortality rate of UK butterflies, possibly due to the increased incidence of diseases and fungi (Radchuk et al., 2013) or decoupling from photoperiod cues (Wiklund et al., 1996). More studies on the wild passerine (a perching bird) in France (Marrot et al., 2017), the great tit (a songbird) in the UK (Cole et al., 2021), and bee colonies across Austria (Switanek et al., 2017) have indicated a significant ecological influence of relatively high temperatures. Hence, from an ecological point of view, a tool that is able to identify CDHW events across all seasons is badly needed.

Zscheischler et al. (2020) proposed a typology of compound hazard events focusing on four themes: preconditioned, multivariate, temporally compounding, and spatially compounding events. For the multivariate events, in which multiple hazards occur at the same time (e.g., droughts and heatwaves in this study), existing literature generally identifies such compound events as the intersection of the periods of the co-occurring hazards (Zscheischler and Seneviratne, 2017; Ridder et al., 2018; Mukherjee and Mishra, 2021). However, from a more systematic point of view, we can distinguish various ways to determine the start and end dates of such compound events: the mentioned intersection, as well as the union or one of the periods conditioned on the other.

The overall objective of this study is to develop a method to identify CDHW events on a daily scale across four seasons by addressing the following research questions:

• a.

  Can we consistently define droughts and heatwaves?

• b.

  How can we define heatwaves across four seasons?

• c.

  How can we avoid the fact that the use of indices calculated on a daily scale could produce a large number of small-scale events and/or mutually dependent events when using fixed-threshold-based identification approaches?

• d.

  How should we interpret the term “compound” and how can we retrieve CDHW events from the droughts and heatwaves themselves?

• e.

  Do CDHW events in non-summer seasons have similar characteristics to those in the summer season?

• f.

  Are the temporal and frequency trends observed consistent with existing studies?

This study is structured as follows. In Sect. 2, we describe the method developed and the data used. We make use of 120 years of daily precipitation and temperature data observed at the Belgian Meteorological Institute in Uccle, Brussels-Capital Region, Belgium. In Sect. 3, we present and discuss the results obtained. Finally, this discussion is followed by some conclusions in Sect. 4.

2 Materials and methods

This section mainly describes the proposed method for the identification of CDHW events, which includes three steps. First, we define droughts and heatwaves and calculate the indices to address research questions (a) and (b). Second, we identify drought and heatwave events from daily indices to address research question (c). Finally, we identify CDHW events from drought and heatwave events to address research question (d).

2.1 Data

Daily minimum and maximum temperatures from 1901 to 2020 were acquired from the climatological station of the Royal Meteorological Institute (RMI) of Belgium in Uccle (50°47^′55^′′ N, 4°21^′29^′′ E; 100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), while daily rainfall was obtained from 10 min precipitation time series recorded at the same site. These high-quality observations have been used in many studies (see, e.g., Verhoest et al., 1997; Ntegeka and Willems, 2008; Vandenberghe et al., 2010; Pham et al., 2018).

2.2 Daily drought and heatwave indices

In this study, a heatwave is defined in a relative sense, i.e., as a period of excessively high temperatures compared with the expected normal temperatures in this period. As the expected normal temperatures differ for individual calendar days, heatwaves defined in this manner could capture anomalous warm periods during summer and winter. We use the daily SPI (McKee et al., 1993) and SHI (standardized heatwave index; Raei et al., 2018) as indices to identify droughts and heatwaves, respectively. The SPI (SHI) is a standardized index to quantify the extent to which precipitation (temperature) deviates from the climatological average. SPI (SHI) describes dry (hot) spells according to the probability of occurrence in a reference period. Given the nonstationarity of temperature records due to global warming, the historical reference period and climate normals for calculating such probabilities must be carefully selected. Accounting for the fact that people and ecosystems adapt to a changing climate, we use the past 30 years as the historical reference period and average the values on every day during this period. This average is then regarded as the expected normal temperature on this day, instead of the average during the longest climatology (the period of record). This 30-year moving-window approach has been suggested to deal with the climate nonstationarity bias (Hoylman et al., 2022). To examine the impact of this approach, we apply the Mann–Kendall (MK) test for the daily mean temperature in the period from 1901 to 2020 (the period of record) and in the 30-year moving windows (1901–1930, 1902–1931,…,1991–2020; 91 windows per day) based on the data described in Sect 2.1. The test results (see Fig. S1 in the Supplement for details) point out that, over the period of record, 40.27 % of calendar days (147 out of 365 d) show a significant increase, while not a single day shows a significant decrease. However, for the 30-year moving windows, only 5.59 % of windows (33215 windows in total, 91 windows per day) show a significant trend, 4.38 % with an increase and 1.21 % with a decrease, indicating that this approach better accounts for climate nonstationarity.

After having overcome the nonstationarity problem in the above way, we calculate SHI values for each day as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\widehat{T}}_{m,i}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{N_{\mathrm{h}}}{\sum }_{j=i-N_{\mathrm{h}}+\mathrm{1}}^i{T}_{m,j},& \text{if }i\ge N_{\mathrm{h}}\\ \frac{\mathrm{1}}{N_{\mathrm{h}}}({\sum }_{j=\mathrm{365}-N_{\mathrm{h}}+i+\mathrm{1}}^{\mathrm{365}}{T}_{m-\mathrm{1},j}& \\ +{\sum }_{j=\mathrm{1}}^i{T}_{m,j}),& \text{if }i<N_{\mathrm{h}};\end{array}\right.\text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{SHI}}_{m,i}={\mathrm{\Phi }}^{-\mathrm{1}}\left(F_{m,i}^{\mathrm{h}}\left({\widehat{T}}_{m,i}\right)\right).\end{array}$$

Here, i and j represent calendar days ($i,j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{365}\mathit{\}}$); $m\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },M\mathit{\}}$ denotes the year considered, where M is the number of years of historical observations; and T_m,i (°C) is the daily mean temperature at day i of year m, which is obtained by averaging the daily maximum and minimum temperatures. For leap years, we average the daily mean temperature on the 28 and the 29 February, assign it to the 28 February, and then remove the leap days. N_h (d) is the length of the accumulation period for a heatwave; ${\widehat{T}}_{m,i}$ (°C) is the mean temperature over the accumulation period up to and including day i; $F_{m,i}^{\mathrm{h}}$ is the cumulative distribution function (CDF) fitted to $\mathit{\{}{\widehat{T}}_{m,i},{\widehat{T}}_{m-\mathrm{1},i},\mathrm{\dots },{\widehat{T}}_{m-\mathrm{29},i}\mathit{\}}$ if m≥30 or to $\mathit{\{}{\widehat{T}}_{\mathrm{30},i},{\widehat{T}}_{\mathrm{29},i},\mathrm{\dots },{\widehat{T}}_{\mathrm{1},i}\mathit{\}}$ otherwise; Φ⁻¹ is the inverse standard normal CDF; and SHI_m,i is the estimated SHI value at day i of year m.

While not relevant for SHI values, a critical problem in the calculation of SPI values is how to deal with the common zeros in precipitation time series (Mishra and Singh, 2010). Indeed, strictly positive distributions (like the gamma distribution) are undefined at zero, while the density function at zero is always zero for a great deal of positive distributions (like the Weibull distribution), no matter how many zero values are present in the data to which they are fitted. If these distributions were used to fit precipitation time series with many zeros, zero and small rainfall amounts (usually indicating dry conditions) would not be well accounted for, which might affect the drought analysis. To resolve this issue, we split the precipitation values into zero values and positive values and estimate their probabilities separately (Naresh Kumar et al., 2009; Farahmand and AghaKouchak, 2015), as expressed in Eq. (4). In this way, we calculate SPI values for each day as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathrm{SP}}_{m,i}=\left\{\begin{array}{ll}{\sum }_{j=i-N_{\mathrm{d}}+\mathrm{1}}^i{P}_{m,j},& \text{if }i\ge N_{\mathrm{d}}\\ {\sum }_{j=\mathrm{365}-N_{\mathrm{d}}+i+\mathrm{1}}^{\mathrm{365}}{P}_{m-\mathrm{1},j}& \\ +{\sum }_{j=\mathrm{1}}^i{P}_{m,j},& \text{if }i<N_{\mathrm{d}};\end{array}\right.\text{(4)}& {\displaystyle }{q}_i={\displaystyle \frac{Z_i}{M}};\text{(5)}& {\displaystyle }{F}_i^{\mathrm{d}}\left({\mathrm{SP}}_{m,i}\right)=\left\{\begin{array}{ll}{q}_i,& \text{if }{\mathrm{SP}}_{m,i}=\mathrm{0}\\ q_i+(\mathrm{1}-q_i)\times F_i\left({\mathrm{SP}}_{m,i}\right),& \text{if }{\mathrm{SP}}_{m,i}>\mathrm{0};\end{array}\right.\text{(6)}& {\displaystyle }{\mathrm{SPI}}_{m,i}={\mathrm{\Phi }}^{-\mathrm{1}}\left(F_i^{\mathrm{d}}\right({\mathrm{SP}}_{m,i}\left)\right).\end{array}$$

Here, P_m,i (mm) is the daily precipitation at day i of year m. For leap years, we add the precipitation on leap days to that on the 28 February and then remove the leap days. N_d (d) is the length of the accumulation period for a drought; SP_m,i (mm) is the precipitation over the accumulation period up to and including day i; Z_i represents the number of zero SP_m,i over the M years; F_i is the CDF fitted to the positive precipitation values $\mathit{\{}{\mathrm{SP}}_{m,i}\mid {\mathrm{SP}}_{m,i}>\mathrm{0},m=\mathrm{1},\mathrm{2},\mathrm{\dots },M\mathit{\}}$; and $F_i^{\mathrm{d}}$ is the CDF of all precipitation values. The SPI value SPI_m,i at day i of year m is then obtained by applying the inverse Φ⁻¹ of the standard normal CDF.

The choice of the distribution to be fitted is an important source of uncertainty in the SPI values (Mishra and Singh, 2010), given that, depending on the case considered, other commonly used distributions are valid (Laimighofer and Laaha, 2022). In this study, we consider 10 commonly used distributions from which we select the option with the lowest value for the Akaike information criterion (AIC) (Akaike, 1974) in the calculation of both the SPI and SHI values. The distributions considered are the normal, exponential, gamma, generalized extreme value (GEV), inverse Gaussian, logistic, log-logistic, log-normal, Burr, and extreme value (EV) distributions.

2.3 Extreme and independent events

Droughts (heatwaves) are commonly identified by a fixed-threshold approach as periods with index values consecutively below (above) a certain threshold (e.g., −1 for SPI), even though some studies use different thresholds for the onset and ending of an extreme event. In this study, we select one threshold and refer to it as the pre-identification threshold and to the resulting events as dry (warm) spells.

However, this fixed-threshold approach, especially when applied to daily indices, may result in a large number of small-scale and mutually dependent spells, causing the identified droughts or heatwaves to be not extreme or not independent. In turn, this could affect the analysis of temporal trends and frequencies. Procedures for excluding minor spells and merging mutually dependent ones have already been explored in the literature (Fleig et al., 2006). However, it is still tricky to discern objectively whether spells are minor or whether two neighboring spells are dependent.

Figure 1 Conceptual illustration of the identification of droughts and heatwaves. Panels (a) and (c) show the dry spells and warm spells obtained after applying the pre-identification threshold; panels (b) and (d) show the final identified droughts and heatwaves after applying the removal and merging procedures. Duration refers to the length of a warm spell or dry spell, while inter-arrival time is the length of the period between the arrivals of two spells. The blue area represents the proximity, while the yellow area represents the severity.

In this study, we develop a method to obtain objective thresholds for excluding minor spells and merging mutually dependent ones. First, two indicators are selected: the duration (d) and the proximity. The duration expresses how long a spell lasts (and equals the number of days in a dry, D_d, or a warm, D_h, spell), while the proximity describes how close two neighboring spells are (Fig. 1). To define the proximity, we first define the total deficiency. The total deficiency of a time interval [a,b] corresponds to the area enclosed between the SPI (or the SHI) curve and the pre-identification threshold line. The total drought deficiency (TD_d) and total heatwave deficiency (TD_h) are calculated as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{\mathrm{TD}}_{\mathrm{d}}=\sum _{i=a}^b({\mathrm{SPI}}_{m,i}-{\mathrm{SPI}}_{\mathrm{d}}),\text{(8)}& {\displaystyle }& {\displaystyle }{\mathrm{TD}}_{\mathrm{h}}=\sum _{i=a}^b({\mathrm{SHI}}_{\mathrm{h}}-{\mathrm{SHI}}_{m,i}).\end{array}$$

Here, SPI_d and SHI_h are the corresponding pre-identification thresholds.

The drought proximity (C_d) or the heatwave proximity (C_h) is defined as the total drought deficiency (TD_d) or heatwave deficiency (TD_h), respectively, of the time window between two neighboring spells. The proximity provides a comprehensive view of the period between two spells, accounting for both duration and how far it is from being dry (or warm). For example, a small value of C_d means that the interval is short and/or the SPI values are near the pre-identification threshold. In that case, the neighboring spells are more likely to be mutually dependent, suggesting that one should merge the neighboring spells and the interval between them into one longer spell.

More importantly, in view of our aim that droughts and heatwaves should be extreme and independent, we assume that their severity should follow a generalized extreme value (GEV) distribution (to guarantee that only extreme-impact events are included) and that their arrival should occur according to a Poisson process (to guarantee that they occur independently in time). Severity is the characteristic used to quantify droughts or heatwaves, usually calculated by means of drought or heatwave indices, respectively. In this study, the severity of a dry (warm) spell is calculated as the negative of the sum of the SPI values (the sum of the SHI values) in a dry (warm) spell. The above hypotheses could be used as objective criteria for determining appropriate thresholds for removing and merging spells, i.e., thresholds that make the severity of the final identified events follow a GEV distribution and concurrently make the inter-arrival times follow an exponential distribution.

Concretely, we remove a spell when its duration is shorter than a chosen removal threshold (D_d<R_d or D_h<R_h). After removing such spells, we check whether any proximity is smaller than a chosen merging threshold (C_d<M_d or C_h<M_h). If this is the case, the two neighboring spells are considered to belong to one longer spell. Therefore, they are merged, and the period between both initial spells is now considered to be part of a longer spell. These processes are illustrated in Fig. 1.

The specific steps to obtain the proper thresholds for removing and merging spells, using the identification of heatwaves as an example, are as follows:

1. Given an accumulation period length N_h and pre-identification threshold SHI_h, we calculate the daily SHI values using Eqs. (1) and (2).

2. We label all warm spells potentially belonging to a heatwave (days consecutively above SHI_h).

3. We consider all possible combinations of removal and merging thresholds R_h and M_h, i.e.,$R_{\mathrm{h}}\in \mathit{\{}\mathrm{1},\mathrm{\dots },\mathrm{5}{N}_{\mathrm{h}}\mathit{\}}$ and $M_{\mathrm{h}}\in \mathit{\{}-\mathrm{\infty },\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{5}{N}_{\mathrm{h}}\mathit{\}}$, where −∞ means that no merging is carried out. For each combination the following substeps are carried out.

   • 3.1

     We remove the warm spells with a duration shorter than R_h and relabel the remaining warm spells. Subsequently, we merge neighboring spells if the corresponding proximity C_h is less than M_h and again relabel the warm spells. After removing and merging, the severities (sum of SHI values during a warm spell) and inter-arrival times (in days) for this combination are calculated as shown in Fig. 1.

   • 3.2

     We test whether the inter-arrival times follow an exponential distribution by means of a Kolmogorov–Smirnov test (KS test) (Massey, 1951) to satisfy the hypothesis that heatwaves are independent. The KS test is chosen because it does not assume any particular underlying distribution. If the KS test rejects the null hypothesis, we proceed to the next combination; otherwise, we go to step (3.3).

   • 3.3

     We test whether the GEV distribution fits the severities to satisfy the hypothesis that heatwaves should be extreme. If the KS test rejects the null hypothesis, we proceed to the next combination; otherwise, we go to step (3.4).

   • 3.4

     We check whether other commonly used distributions also fit the severities by means of the KS test. The distributions considered are the same as those used in the calculation of the SPI and SHI values. If the KS test rejects all other distributions, we make a note that the GEV distribution is the only fit for this combination and then proceed to the next combination; otherwise, we go to step (3.5);

   • 3.5

     We compute the AIC and root-mean-square deviation (RMSD) values for all distributions that pass the KS test. We record whether the GEV distribution has the lowest AIC or the lowest RMSD.

4. Among all combinations for which the GEV distribution is the only fit or results in both the lowest AIC and the lowest RMSD values, we select the one leading to the largest number of events in order to increase the probability of having more compound events and, thus, benefit further analyses. If there is no such combination, we consider the combinations resulting in the lowest AIC only. If no combination at all has been identified, the accumulation period length and pre-identification threshold considered are discarded.

5. If a combination of thresholds has been identified, the warm spells obtained after applying the removal and merging procedures are finally labeled as heatwaves, satisfying the hypothesis that they are extreme and independent.

The identification of drought events follows a similar procedure.

Figure 2 Conceptual diagram of four types of CDHW events: intersection (d-and-h), union (d-or-h), drought conditioned on heatwave (d-cond-h), and heatwave conditioned on drought (h-cond-d). The gray hatching represents the days in droughts or heatwaves; arrows indicate the duration of CDHW events.

2.4 Identification of compound events

After having identified droughts and heatwaves, we proceed with the identification of CDHW events. In the literature, a CDHW event is generally defined as the period of overlap, i.e., the intersection, of a drought and a heatwave (Fig. 2d). However, this is not the only way in which a compound event could be defined. In general, we can distinguish at least four different ways to define a compound event. These compound events consist of the following:

• i.

  Days belonging to both the drought and the heatwave, i.e., the intersection (d-and-h). Note that a drought might overlap with multiple heatwaves, resulting in multiple compound events, and the same applies to a heatwave.

• ii.

  Days belonging to the drought or the heatwave, i.e., the union (d-or-h), provided the events overlap. Note that a drought might overlap with multiple heatwaves, resulting in multiple overlapping compound events that need to be subsequently merged.

• iii.

  Days belonging to the drought, denoted as d-cond-h, provided the events overlap. Note that a drought overlaps with multiple heatwaves. Thus, in the d-cond-h case, only those droughts are preserved that overlap with at least one heatwave.

• iv.

  Days belonging to the heatwave, denoted as h-cond-d, provided the events overlap. Note that a heatwave might overlap with multiple droughts. A heatwave that does not overlap with at least one heatwave is not considered in the h-cond-d case.

Figure 2 displays a conceptual diagram of the above four different types of CDHW events. A CDHW event is characterized by three variables: the duration, the marginal drought severity, and the marginal heatwave severity. The marginal drought severity is computed as the negative of the sum of the SPI values within the CDHW event, while the marginal heatwave severity is computed as the sum of the SHI values within the CDHW event.

It could be of relevance to ecologists to investigate whether compound events have a higher impact than droughts or heatwaves only. For certain species, it is known that heatwaves have an impact (Pipoly et al., 2022), but it might be interesting to investigate the possible intensified or weakened impacts of heatwaves during drought conditions (i.e., h-cond-d). Similarly, droughts occurring during a heatwave may have an increased impact (i.e., d-cond-h). For species that suffer during both droughts and heatwaves, a compound event during which a drought and/or a heatwave occur (i.e., d-or-h) may have a stronger impact than when droughts or heatwaves occur independently of each other.

2.5 Probabilities of occurrence

The probabilities of occurrence of a day being in a drought (P_d), a heatwave (P_h), or a CDHW event (i.e., P_d-and-h, P_d-or-h, P_d-cond-h, and P_h-cond-d) are estimated empirically as the fractions of days belonging to the corresponding events within the study period. By comparing P_d-and-h with the probability of having a compound event assuming that droughts and heatwaves occur independently of each other, i.e., by comparing P_d-and-h with P_dP_h, we can examine whether there is a positive dependence (P_d-and-h>P_dP_h ), negative dependence (P_d-and-h<P_dP_h), or no dependence at all (P_d-and-h≈P_dP_h) between droughts and heatwaves (Ridder et al., 2020). Further, we calculate the likelihood multiplication factor (LMF; Ridder et al., 2020) as $\mathrm{LMF}=\frac{P_{\text{d-and-h}}}{P_{\mathrm{d}}\times P_{\mathrm{h}}}$. In the case of independence, the LMF approximately equals 1; if droughts and heatwaves are positively dependent, the LMF exceeds 1 and increases with the strength of dependence.

Figure 3 The frequency of the distribution selected with the lowest AIC value in the calculation of the SPI and SHI.

3 Results and discussion

3.1 SPI and SHI

We calculated the SPI and SHI for various accumulation period lengths (N_d = 15, 30, 45, 60, 90 d; N_h = 3, 5, 7, 10, 15 d). The distributions in the SPI and SHI are selected based on the lowest AIC values. The frequency of the selected distribution is shown in Fig. 3. The gamma distribution is the most frequently selected one in the calculation of SPI, followed by the GEV, Burr, and normal distributions. For SHI, there is a wider range of distributions, including the normal, gamma, GEV, inverse Gaussian, EV, logistic, log-logistic, and Burr distributions.

The choice of distribution can significantly influence the SPI and SHI values, especially in the tails of the distribution, which are particularly relevant for extreme events. For example, we conducted a comparison between values calculated using fixed distributions (commonly chosen: gamma and GEV for SPI and normal and GEV for SHI) and the values in our study, depicted in Fig. S2. We observed several problematic SHI values when using the GEV distribution, as indicated by the gray circle. Existing studies, such as Mishra and Singh (2010) and Laimighofer and Laaha (2022), have also highlighted the impact of the choice of the distribution to be fitted. We suggest choosing the proper distribution according to the goodness-of-fit criterion.

Table 1 Removal thresholds (R_d), merging thresholds (M_d), average number of days per year in an event, and average number of events per year for (a) droughts and (b) heatwaves. The italic numbers indicate the SPI_d or SHI_h values.

3.2 Removal and merging thresholds

We identify the droughts and heatwaves for various accumulation period lengths (N_d = 15, 30, 45, 60, 90 d; N_h = 3, 5, 7, 10, 15 d) and pre-identification thresholds (SPI_d = −0.5, −1, −1.3; SHI_h = 0.5, 1, 1.3). Table 1 shows the resulting removal and merging thresholds, the corresponding average number of days per year in an event, and the average number of events per year. We consider, for example, the case N_d = 15 d and SPI_d = −1. In this case, the combination R_d = 14 d and M_d = 15, results in 28.3 drought days per year and 1.13 drought events per year on average. The values R_d = 14 d and M_d = 15 indicate that all dry spells shorter than 14 d are removed and that all neighboring spells with C_d≤15 are merged, respectively, to satisfy the hypothesis that drought events are extreme and independent.

The proposed removal and merging method shows the desired behavior. As could be expected, more extended accumulation periods (i.e., larger N_d or N_h values) cause the SPI and SHI values to evolve more gradually over time, resulting in dry or warm spells being longer. Moreover, larger pre-identification thresholds for droughts (i.e., larger SPI_d) or smaller pre-identification thresholds for heatwaves (i.e., smaller SHI_h) tend to lead to the inclusion of more days in dry or warm spells, hence also resulting in their duration being longer. We can expect the removal threshold needed to exclude minor spells to increase with larger N_d or N_h values as well as with larger SPI_d or smaller SHI_h values.

As shown in Table 1, for example, if SPI_d = −1, then R_d takes values of 14, 22, 31, 45, and 59 d when N_d ranges over 15, 30, 45, 60, and 90 d. If N_d = 15, then dry spells shorter than 10, 14, and 23 d are removed when SPI_d ranges over −1.3, −1, and −0.5. In the same way, if SHI_h = 1, then R_h takes values of 9, 10, 11, 13, and 17 d when N_h ranges over 3, 5, 7, 10, and 15 d. In the case that N_h = 3, warm spells shorter than 6, 9, and 13 d are removed when SHI_h ranges over 1.3, 1, and 0.5. The same behavior is also observed for other accumulation period lengths and other pre-identification thresholds. In summary, the way that R_d and R_h change with varying accumulation period lengths and pre-identification thresholds in this study is consistent with expectations.

In contrast to the removal thresholds, no clear patterns can be distinguished for the merging thresholds when varying the accumulation period lengths or pre-identification thresholds in Table 1. Interestingly, the merging thresholds influence the number of events. What stands out in Table 1 is the case in which SHI_h = 1.3 and N_h = 15 d. It then holds that M_h = 11, leading to 0.58 heatwaves per year on average, which is considerably less than in other cases with the same accumulation period length or pre-identification threshold. This can be explained by the fact that a larger merging threshold results in neighboring spells being combined, ultimately leading to fewer events.

Another interesting observation in Table 1 is that no merging is carried out in some cases (M_d = −∞ or M_h = −∞). This is because, using the proposed method, the removal procedure is carried out first, followed by the merging procedure. After the minor spells are removed, the inter-arrival time between spells correspondingly becomes larger, and the neighboring spells are less likely to be mutually dependent. Therefore, it could happen that, after removing minor spells, the inter-arrival time already follows an exponential distribution. Thus, there is no further need to carry out the merging procedure.

We also apply the fixed-threshold method to see what would happen if we did not carry out the removal and merging procedures. The numbers of days and events are strongly reduced when using our method compared with the fixed-threshold one, as only extreme and independent events are selected. We find an average 28.3–33.6 and 17.9–24.9 d yr⁻¹ in droughts for SPI_d = −1 and SPI_d = −1.3, respectively. Without the removal and merging procedure, there are 56.8–60.6 and 35.0–38.4 d yr⁻¹ for SPI_d = −1 and SPI_d = −1.3, respectively (see Table S1 in the Supplement for details). Similarly, there are 19.2–33.6 d yr⁻¹ in heatwaves when SHI_h = 1, whereas there are 68.2–68.8 d yr⁻¹ in heatwaves before removing and merging. This way, only about half of the days are left in extreme events. Moreover, if the fixed-threshold method is applied, 35.5 %–50.7 % of dry spells last less than 5 d and 5.6 %–22.3 % of warm spells only last 1 d in all cases (Fig. S3). We doubt that such short spells influence society or ecosystems and, thus, suggest excluding them from the analysis of climate extremes.

3.3 Identification results for CDHW events

Different accumulation period lengths and/or pre-identification thresholds for droughts and heatwaves could result in different CDHW events. We consider, for instance, the CDHW events corresponding to N_d = 15 d and SPI_d = −1 and to N_h = 3 d and SHI_h = 1. The corresponding removal and merging thresholds are given by R_d = 14 d and M_d = 15 and by R_h = 9 d and M_h = 0.

Next, we discuss the results for the four different types of compound events. In Uccle, in the period from 1901 to 2020 (the period of record), 33 CDHW events of types d-and-h and h-cond-d are identified in total, while 29 CDHW events of types d-cond-h and d-or-h are identified. Hence, their frequency of occurrence ranges from once every 3.6 years to once every 4.1 years (3.6- to 4.1-year return period). Furthermore, the corresponding numbers of days in these CDHW events are given by 2.9 d yr⁻¹ (type d-and-h), 9.1 d yr⁻¹ (type d-or-h), 8.3 d yr⁻¹ (type d-cond-h), and 3.7 d yr⁻¹ (type h-cond-d) on average. The low numbers of events and numbers of days illustrate that the identification of CDHW events only considers rare events.

Figure 4 Identification results in 1976 for droughts, heatwaves, and CDHW events (of types d-and-h, d-or-h, d-cond-h, and h-cond-d) based on the daily SPI (left) and 1-month SPI (right). Gray areas indicate CDHW events. Identification results in other years are shown in Fig. S4.

Given that we work on a daily scale, a more precise identification of the start and end dates of drought events is possible, prompting a more precise identification of CDHW events. To demonstrate this, we take the year 1976 as an example, in which extraordinary droughts and heatwaves occurred in Belgium, leading to several forest fires and even requiring restrictions on the use of potable water (KMI, 2017). Our method identifies a heatwave spanning from 23 June to 10 July, lasting more than half a month, which is exceptionally long for a heatwave in Belgium (shown on the left of Fig. 4). At the same time, there was a severe drought event lasting more than 4 months, starting on February 28 and ending on 11 July (due to an intense precipitation event – 31.7 mm daily rainfall – on 12 July). Moreover, a CDHW event of types d-or-h and d-cond-h is identified starting on 28 February and ending on 11 July, while a CDHW event of the other two types is identified from 23 June to 10 July. Applying a 1-month SPI (shown on the right of Fig. 4), the drought is identified as ending on 30 June. Consequently, the corresponding compound event of type d-and-h then starts on 23 June and ends on 30 June, lasting only 8 d, in contrast to 18 d if the daily SPI is used. Comparison of these two results reveals that the use of the monthly index could significantly bias the estimation of the duration of CDHW events. The proposed method offers a more precise identification, as it is able to capture the start and end dates of droughts and CDHW events on an intra-monthly timescale thanks to the use of daily indices.

Moreover, the method proposed also allows one to identify heatwaves and CDHW events in the winter season because the SPI and SHI treat each day in the same way by evaluating extreme anomalies. Figure 5 shows an example of the winter of 2015 when Uccle had an exceptionally mild winter. An average temperature of 6.6 °C was recorded from November 1 until the end of December, while the average temperature in this period over the past 30 years was 3.8 °C. A CDHW event of type d-and-h is identified from 3 to 17 November. Given that CDHW events in non-summer months are also of interest to ecologists (Williams et al., 2010), our method could be of high interest to them.

Figure 5 Identification results for a winter event in 2015. The orange color is used for droughts, while the green color is used for heatwaves.

Figure 6 Scatterplots for the duration and marginal severity of CDHW events of four types. r stands for the Pearson correlation coefficient, with the subscripts 1, 2, and 3 representing the duration, marginal drought severity, and marginal heatwave severity, respectively. The orange color is used for droughts, while the green color is used for heatwaves.

CDHW events can be described by the following numerical characteristics: the duration, the marginal drought severity, and the marginal heatwave severity. The pairwise relationships between these three characteristics are explored (see Fig. 6): all show positive correlations for the four types of CDHW events. Especially for type d-and-h, the three characteristics are highly correlated, with Pearson correlation coefficients larger than 0.9. Interestingly, CDHW events of type d-or-h and d-cond-h tend to have a longer duration (with a maximum duration of 148 d), while the duration of CDHW events of type h-cond-d and d-and-h is much shorter (with a maximum duration of 32 and 28 d, respectively), because the length of heatwaves limits the duration of these events.

3.4 Temporal trends

In this subsection, we specifically analyze the CDHW events corresponding to N_d = 15 d and SPI_d = −1 and to N_h = 3 d and SHI_h = 1. However, the trends observed are in line with those for other accumulation period lengths and pre-identification thresholds (see Table S2).

Figure 7 Temporal trends in droughts, heatwaves, and CDHW events of four types. Two periods are considered: 1901–1960 (gray) and 1961–2020 (yellow). The superimposed text indicates the number of days in these events per year and the number of events per year averaged over the corresponding period.

We can observe that the number of CDHW events and the number of days in CDHW events in Uccle are increasing over the period from 1901 to 2020 (see Fig. 7). In particular, the number of CDHW events shows a large increase from 0.167 events per year in the period from 1901 to 1960 to 0.317 events per year in the period from 1961 to 2020 for types d-or-h and d-cond-h (or from 0.217 to 0.333 events per year for types d-and-h and h-cond-d), which is nearly double. Also the numbers of days in the four types of CDHW events are higher in the period from 1960 to 2020 compared with the period from 1901 to 1960, such as a rise from 7.78 to 10.4 d yr⁻¹ for CDHW events of type d-or-h and from 2.23 to 3.50 d yr⁻¹ for CDHW events of type d-and-h. We further applied the MK test to check whether there is a statistically significant trend in the number of days per year in CDHW events. Despite the positive trends for the CDHW events of all four types considered, they do not appear to be significant at the 95 % confidence level.

The increasing frequency of heatwaves contributes to the increasing frequency of CDHW events. It should be underlined that heatwaves are obtained through comparison with the expected normal temperatures during the past 30 years instead of the period of record. Nevertheless, there has been a substantial rise in the number of days in heatwaves per year and the number of heatwaves, increasing from 15.6 d yr⁻¹ and 76 events in the period from 1901 to 1960 to 22.8 d yr⁻¹ and 111 events in the period from 1961 to 2020. In addition, the MK test for the number of days in heatwaves per year confirms a significant increasing trend at the 95 % confidence level. For droughts, there is a slight but nonsignificant increase in the number of days in droughts per year, despite the rise in the average number of days in droughts per year and the number of droughts from 27.3 d yr⁻¹ and 63 events in the period from 1901 to 1960 to 29.3 d yr⁻¹ and 72 droughts in the period from 1961 to 2020. Thus, the increases for the CDHW events are mainly due to the corresponding increases for heatwaves. This is consistent with the study of Manning et al. (2019), who analyzed dry and hot events over Europe in the summer season.

Figure 8 Seasonal characteristics. The first three rows show the probabilities of occurrence of a day being in a drought, a heatwave, or a CDHW event of types d-and-h, d-or-h, d-cond-h, and h-cond-d in each month and each season. The bottom two panels compare P_d-and-h to P_dP_h in each month and each season, respectively.

3.5 Seasonality

Next, we examine the seasonality of CDHW events (Fig. 8). Their probabilities of occurrence show monthly characteristics. For instance, CDHW events in April, August, and October all have relatively high probabilities of occurrence compared with the other months, regardless of the type of CDHW event considered (d-and-h, d-or-h, d-cond-h, or h-cond-d). Contrarily, the probability of experiencing a compound event is low in December, January, and February. For the month of January in particular, no CDHW events have been identified in the period from 1901 to 2020. The limited number of events may preclude the reliable estimation of the monthly probability of occurrence. Alternatively, we can also calculate the probability of occurrence on a seasonal scale. The results show that the probability of occurrence of CDHW events in winter is smaller than in the other three seasons, but this difference is not apparent for droughts or heatwaves.

We further compare the probability of occurrence of CDHW events of type d-and-h (P_d-and-h) and the probability of having a compound event assuming that droughts and heatwaves are independent (P_dP_h). As can be seen from Fig. 8, P_d-and-h is smaller than P_dP_h in winter, whereas it is greater in spring, summer, and autumn. This reveals that droughts and heatwaves in Uccle tend to coincide in spring, summer, and autumn, while CDHW events are less likely in winter.

Atmospheric and land–atmosphere interactions could explain the above seasonal differences. It has been widely reported and interpreted in the literature (see, e.g., Geirinhas et al., 2021; Schumacher et al., 2019; Miralles et al., 2014) that droughts and heatwaves are positively dependent in summer. First, this phenomenon has been interpreted in the context of soil moisture–atmosphere coupling, with soil moisture deficits related to reduced precipitation leading to enhanced surface sensible heating and higher surface temperatures (Berg et al., 2015). Second, heatwaves could reduce the total energy transfer to the atmosphere, possibly decreasing convective precipitation (Zaitchik et al., 2006). This, in turn, leads to a soil–precipitation feedback loop that tends to extend or intensify drought conditions (Mazdiyasni and AghaKouchak, 2015). Third, the synoptic-scale weather systems favorable for extreme heat are also unfavorable for rain, which increases the probability of co-occurrence (Berg et al., 2015). Precipitation in Uccle is characterized by cyclonal and convective rainfall in the summer months (Verhoest et al., 1997). The above three mechanisms are possible explanations for the simultaneous occurrence of CDHW events in the non-winter seasons at Uccle. In winter, negative correlations between droughts and heatwaves dominate in Uccle, as warm moist advection in extratropical cyclones favors precipitation (Trenberth and Shea, 2005). Uccle is located in the westerlies area, with the prevailing wind direction from west to east. The westerlies are strongest in the winter hemisphere; the pressure is lower over the poles during this time. The warm winter west wind in Uccle favors precipitation, causing CDHW events to be less frequent. Furthermore, due to the radiative influence of clouds, cloudy (often wet) days in winter are generally warmer, whereas sunny (dry) days are colder. Finally, the development of CDHW events is less likely in winter, as soil–precipitation feedbacks play no role because the low atmospheric demand can easily be met by the typically wet soils (Seneviratne et al., 2010; Schumacher et al., 2022).

3.6 A small validation experiment

As is well known, it is challenging to validate the identification of drought, heatwave, and CDHW events due to the absence of a universally accepted definition as well as the difficulty in measuring the exact start and end dates. We designed an experiment to validate the identification indirectly and compare our proposed method with the commonly used method in other studies. The experiment is based on the well-established physical mechanism that a strong positive land–atmosphere feedback exists between droughts and heatwaves during summer. Therefore, droughts and heatwaves are highly dependent in summer, generating a higher probability of CDHW events than when droughts and heatwaves are independent.

The identification method usually comprises two (or three) steps: pre-identification and removal (with merging in some studies). These steps categorize spells into two distinct groups:

• A:

  final identified drought and heatwave events;

• B:

  removed dry and warm spells.

Furthermore, when coupled with known physical mechanism, if the identification method is effective, we can expect that, during summer, the removed spells in group B would exhibit less dependence than the events in group A. Employing the LMF to quantify the strength of the dependence, we expect the following behavior during summer:

1. LMF_A should be greater than 1;

2. LMF_A should be greater than LMF_B.

We consider all possible combinations of different accumulation period lengths and/or pre-identification thresholds for droughts and heatwaves – N_d = 15, 30, 45, 60, and 90 d; N_h = 3, 5, 7, 10, and 15 d; SPI_d = −0.5, −1, and −1.3; and SHI_h = 0.5, 1, and 1.3 – resulting in 5 × 5 × 3 × 3 = 225 scenarios in total. For each scenario, we calculated LMF_A, LMF_B, and LMF_A−LMF_B in June, July, and August for our proposed removal and merging method as well as for the commonly used method in existing studies by employing a subjective and fixed minimum duration to exclude minor spells. For example, many studies retain drought events only if they persist for more than 30 d to focus on significant events (Brunner and Stahl, 2023; Xu et al., 2023; Christian et al., 2021). Similarly, heatwave events are commonly identified as lasting at least 3 consecutive days (Ridder et al., 2020; Yin et al., 2023).

Validation results (Figs. S5, S6) show that LMF_A consistently exceeds 1 for both methods, which follows the aforementioned expected behavior (1). However, a notable distinction arises in expected behavior (2). In all 225 scenarios examined, LMF_A−LMF_B consistently takes positive values, indicating our method's robustness. In contrast, for the fixed removal method, the condition ${\mathrm{LMF}}_{\mathrm{A}}-{\mathrm{LMF}}_{\mathrm{B}}>\mathrm{0}$ is not fulfilled in 36 out of 225 scenarios, i.e., in 16 % of cases, suggesting that the removed dry and warm spell events had a higher degree of dependence than the final identified events.

Overall, applying a fixed removal threshold method for various accumulation periods and pre-identification thresholds of droughts and heatwaves introduces the potential risk of generating unreasonable results. In contrast, our removal and merging method effectively addresses this challenge during validation.

4 Conclusion

We proposed a method to identify droughts, heatwaves, and compound events. The identification on a daily scale systematically and objectively removed minor spells and merged mutually dependent ones. The analysis conducted at Uccle demonstrates the effectiveness of the proposed method in four ways. First, the values of removal thresholds exhibit the desired behavior, adapting effectively to varying accumulation periods and pre-identification thresholds. Second, the frequency of occurrence of heatwaves and CDHW events has increased in the period from 1961 to 2020 compared with the period from 1901 to 1960, and the increasing temperatures contribute to the increase in CDHW events, which aligns with the current literature (Manning et al., 2019). Moreover, the occurrence of CDHW events shows seasonal patterns, with the occurrences of droughts and heatwaves being negatively dependent in winter but positively dependent in the other three seasons, which could be explained by atmospheric and land–atmosphere interactions. Fourth, a validation experiment based on this positive dependence in summer demonstrated the robustness of the proposed method compared to commonly used methods. We used one station to demonstrate the method, but further studies are needed to validate whether the proposed identification method works well in different climatic zones.

By upscaling the temporal resolution to a finer one, our daily-scale identification captures variations that monthly scales often miss, providing more precise event start and end dates for droughts and CDHW events. This more precise identification could enhance the capacity for detection, assessment, monitoring, and early warning of both drought events and CDHW events.

Furthermore, our definition in relative terms allows for the identification of heatwaves and CDHW events across all four seasons, including non-summer periods. This expanded understanding is crucial, as it sheds light on the ecological repercussions that extend beyond the confines of the traditional summer-focused perspective. The ecological impacts of CDHW events in non-summer seasons are also significant. For instance, in regions characterized by temperate continental and temperate monsoon climates, CDHW events in non-summer seasons link to wildfire weather (Tian et al., 2011). In such regions, the winter season itself often represents the dry season, characterized by reduced precipitation and frequent strong winds. The dry season becomes even drier when drought conditions co-occur with abnormally high temperatures. This exacerbates the dryness of the soil and the atmosphere, accelerating the drying of forest litter and setting the stage for an elevated risk of wildfires.