European heatwave tracks: using causal discovery to detect recurring pathways in a single-regional climate model large ensemble

Observational records provide a limited amount of heatwave events in Europe. Thus, employing a SMILE allows to increase the sample size: By using the Canadian Regional Climate Model, version 5, Large Ensemble (CRCM5-LE, Leduc et al 2019) we obtain 50 members of comparable climate statistics to derive robust results at high spatial resolution (0.11^∘, 12.5 km) while considering internal climate variability. While the CRCM5-LE was already extensively used for analyses on rare (hydrometeorological) extremes (Champagne et al 2020, Wood and Ludwig 2020, Böhnisch et al 2021, Brunner et al 2021, Poschlod and Ludwig 2021), the present study is the first to apply causal inference where large samples are also of importance (Spirtes and Glymour 1991, Runge et al 2019b).

We pool the period 1981–2010 from 50 members into one 1500 year time series. Since all members are designed to share comparable climate statistics, typical patterns and subsequently pathways of phenomena are assumed to be comparable across all members (and to reality), only varying due to internal climate variability (Leduc et al 2019).

For heatwave definition, we use daily maximum temperature (tasmax). Daily 500 hPa geopotential height (Z500) is employed for large-scale atmospheric pattern definition. We first perform heatwave clustering and track derivation on a ERA-Interim driven (Dee et al 2011) run of CRCM5 (CRCM5/ERA, Leduc et al 2019) in order to evaluate them against historical events while using the same regional climate model. Additionally, an external evaluation on the E-OBS gridded dataset (version 22.0e at 0.1^∘ spatial resolution; Haylock et al 2008, Cornes et al 2018) is performed (figure S1). All analyses are confined to land areas within a European domain (Leduc et al 2019).

Heatwave definitions in literature vary according to study perspectives (Perkins and Alexander 2013, Smith et al 2013, Shafiei Shiva et al 2019). To properly assess heatwave tracks, our heatwave definition requires information for each time step (day) under consideration. We define heatwaves relative to local climatology (Della-Marta et al 2007) as a minimum of three consecutive hot days, separated by at least three non-hot days (Keellings and Moradkhani 2020, Spensberger et al 2020) from the preceding/following hot day: a hot day occurs if the centered three-day-running mean of tasmax exceeds the local 95th JJA (1981–2010) percentile (similar to Lhotka and Kyselý 2015a). Negative anomalies are omitted. Whereas the percentile in a single run or observations/reanalysis may be poorly estimated due to the limited sample size (30 years), uncertainty due to internal variability is reduced in the SMILE by deriving the percentile from a 1500 year sample. Additionally, percentile based anomalies reduce potential inconsistencies among the data sources induced by climate model bias.

Spurious artifacts ($\leqslant\kern-0.5pt{9} \times 9$ pixels) and isolated hot day occurrences are removed. A time step is considered as belonging to an extended heatwave event, if hot days occur in at least $1\%$ of the land area in the domain (i.e. 500 grid cells). All matching time steps are next used to derive core regions by clustering (Felsche et al 2022). Preliminary analyses revealed hot days in months May–October, such that we splice this period over all years and members (Ebert-Uphoff and Deng 2012a).

For core region definition, we refer to spatially coherent recurring patterns of European heatwaves by means of clustering (Felsche et al 2022).

Clustering is frequently employed for mid-latitude weather pattern analysis and classification (Smyth et al 1997, Stefanon et al 2012, Lhotka and Kyselý 2015a, Hannachi et al 2017, Wang et al 2018, Keellings and Moradkhani 2020, Machado and Lopes 2020). We use daily maps of hot day occurrences as input to a two-step agglomerative hierarchical clustering (Felsche et al 2022). The optimal number of clusters is determined by computing the distortion score for every possible number of clusters and piking the knee of the curve (Jung et al 2003). When comparing cross-validation clustering results against cluster results from a Monte–Carlo pseudo-experiment, nine spatial heatwave patterns are obtained which are significant on the 99% level based on a two-sided t-test (Felsche et al 2022). These cluster footprints correspond well with historical events captured by other heatwave indexes (e.g. Russo et al 2015). Felsche et al (2022) shows high consistency among clusters in CRCM5/ERA and CRCM5-LE. We follow their naming conventions for CRCM5-LE: BI (Britain and Ireland), IP (Iberian Peninsula), WE (Western Europe), CEE (Central-Eastern Europe), SCA (Scandinavia), NEE (North-Eastern Europe), SEE (South-Eastern Europe), GSI (Greece/Southern Italy). Sequences of time steps labeled according to their respective clusters then represent single heatwaves (figure 1).

Zoom In Zoom Out Reset image size

Next, clustered time steps are spatially averaged, using the normalized intensity of temporally averaged cluster footprints as weights for their centroid coordinates (Herrera-Estrada et al 2017, Herrera-Estrada and Diffenbaugh 2020, Sánchez-Benítez et al 2020). This provides nine time series of positive temperature anomalies, located at the cluster footprint centers (i.e. their centroids) and representing the components (Runge et al 2015) among which spatio-temporal links are to be established. A sensitivity experiment using different weights for aggregation and therefore time series positioning revealed no major changes in link directions. Since CRCM5/ERA and CRCM5-LE provide slightly different cluster footprints, the positions of these time series diverge as well (see figure 2(a) and figure S1 for E-OBS).

Zoom In Zoom Out Reset image size

In order to analyze heatwave tracks between heatwave onset and offset, we split them into transitions among core regions. This captures direct movements between various regions. Location changes within cluster regions or growth and shrinkage of the spatial extent without heatwave core shifts are considered as stationary.

A heatwave moves from region A to region B, if heatwave onset in region A is followed by heatwave onset in region B, and if heatwave onset in region B was preceded by heatwave onset in region A: The heatwave stays in region A until day t and moves to region B at day t+1. We thus aim to explain a later position of a heatwave by its previous position (i.e. effect preceded by its cause).

Generally in causal discovery, distinct nodes representing single processes (variables or locations of the same variable) are derived prior to link establishing, by e.g. sampling of equally spaced grid points or rotated principal component analysis (Ebert-Uphoff and Deng 2012b, Deng and Ebert-Uphoff 2014, Runge et al 2015, 2019a, Nowack et al 2020). Here, we employ time series assigned to nine spatially distributed clusters of heatwave core regions.

To define linear spatio-temporal links among heatwave core regions we employ a causal discovery algorithm provided with the Tigramite 5.1 Python package, namely the PCMCI algorithm (Runge et al 2019a, Nowack et al 2020). It is based on the PC algorithm introduced by Spirtes and Glymour (1991): First, a fully connected graph of all variables (nodes) under consideration, e.g. time series at different locations, is created. Then, conditionally independent links (edges) are deleted, while for the remaining links directions are identified (e.g. by considering a temporal ordering of cause and effect) (Ebert-Uphoff and Deng 2012a, Deng and Ebert-Uphoff 2014, Runge et al 2019b). An amelioration with respect to time series is provided by the two-step PCMCI algorithm, adding a momentary conditional independence test to the PC algorithm (Runge et al 2019a, Nowack et al 2020): it rejects spurious or indirect links by testing for independence conditional on the common past of network nodes. Thus, only direct links among network nodes are kept which goes beyond pure pairwise correlation graphs (Runge et al 2019b). The PCMCI algorithm can be more powerful than correlation analysis in determining relationships between variables, even in finding links which may not be obvious in classic correlation analysis (Almendra-Martín et al 2022). Additionally, it bears the potential to find indirect links and consider common drivers as opposed to other forms of causal discovery like Granger causality (Granger 1969, Runge et al 2019b, Galytska et al 2022).

In this present study, conditional independence tests are conducted based on partial linear correlations (ParCorr) with a significance threshold of $\alpha = 10^{-3}$. We use a minimum temporal lag of one day, reflecting transitions, and a maximum of three days, reflecting the three days separating two heatwaves. Additionally, we focus on positively correlated links because negative correlations would correspond to the non-occurrence of heatwaves in region B after occurrence in region A.

The importance of heatwave drivers varies regionally (Zschenderlein et al 2019), with dynamical mechanisms generally dominating over local thermodynamical ones (Suarez-Gutierrez et al 2020). As mentioned previously, high pressure, i.e. positive Z500 height anomalies, may lead to subsidence, cloud dispersal and subsequently enhanced heating at the surface. Additionally, prevailing dry conditions, e.g. low soil moisture contents, amplify the heating process by reducing the latent heat flux. Both drivers may interact and intensify, e.g. in feedback loops (Fischer et al 2007), since prevailing heat in turn may increase evaporation, thereby leading to further soil drying (Schumacher et al 2019).

We focus on associated large-scale Z500 conditions to investigate why heatwaves follow a given track. Patterns similar to the Z500 anomaly composites alongside selected heatwave transitions are searched during months May–October. Thereupon we calculate the probabilities of (a) having this pattern during an eight-day period centered at the indicated transition and (b) of having this transition during the indicated pattern. Pattern similarity is assessed by means of masking: if at least 66% of grid cells at an arbitrary time step show positive anomalies inside the composite positive anomaly region, and if at the same time an analogous agreement is achieved within the composite negative anomaly region, this time step is considered as being similar to the composite pattern. Consecutive days of similar patterns are summarized as one occurrence, which is related to atmospheric circulation persistence (see Jézéquel et al 2018, who relate atmospheric patterns to hot days by means of Euclidean distance calculation).

During 1981–2010, 115 European heatwaves with at least three days are identified in CRCM5/ERA compared to 5425 in the 50-member CRCM5-LE. Most heatwaves occur during summer months JJA, starting earlier in the West than in the East (Felsche et al 2022). The relative frequency of heatwaves ending in the Eastern parts of Europe is increased compared to the frequency of heatwaves starting there (figure 2(a)). Accordingly, heatwaves start more frequently in the Western parts of Europe. Taken together, these first findings indicate a general west-to-east movement of heatwaves in both CRCM5/ERA and CRCM5-LE, which mirrors the dominating westerly flow onto Europe (e.g. Zschenderlein et al 2019). Heatwaves starting outside the defined core regions, e.g. outside the domain, over Northern Africa or the ocean, are not captured.

Most heatwaves consist of several transitions (figures 2(b) and (c)). While in general the number of transitions increases with heatwave length, very stationary heatwaves extend up to 20 days within one cluster and highly mobile heatwaves encounter 6 transitions in ten days. CRCM5-LE shows a stronger variability among transition–length combinations than CRCM5/ERA, representing natural variability of heatwave characteristics. The three historical examples 2003, 2006, and 2010 exhibit extreme lengths with respect to the majority of events (e.g. three transitions in 44 days in 2010).

A more detailed network of statistical robust sequences among neighboring core regions is provided in figures 3(a) and (b). It represents a mixture of latitudinal and longitudinal transition directions like a spine of the continental landmasses: While in the West, a South-Northward direction prevails, the East experiences the opposite direction. Within the scope of our study, we focus on links with lag $= {1}$ among two core regions at a time. In order to obtain robust results, links were first derived on the entire 1500-summer period. To include internal variability of these, they were also calculated for each member separately, allowing to count their frequency. In figure 3, only links occurring in at least 15 members (i.e. 30% of the ensemble) are presented with colors indicating their frequency (compare Galytska et al 2022, Karmouche et al 2022, for alternative presentation methods in ensembles). As a real-world example, the July 2007 heatwave (figure 1) is represented within the SEE–GSI link.

Zoom In Zoom Out Reset image size

In order to further evaluate the obtained links, we next compare the graph based on model data with reanalysis and observation based links (compare Karmouche et al 2022). While most of these links figure in the CRCM5/ERA and E-OBS network as well (figures S1(a) and (c)), several links from historical data are missing within the robust CRCM5-LE links, being thus (most likely) rare outliers. Others are shifted compared to CRMC5-LE since the position of core regions diverges. The CRCM5-LE links show higher correspondence to the representation in CRCM5/ERA (F1-scores among members and reference above 0.5, figure S2) than to the one in E-OBS (F1-scores close to 0.2; most likely due to clearly shifted or split core regions like, e.g. GSI). Additionally, no link connecting CEE and SCA can be established in CRCM5/ERA as was also found in Spensberger et al 2020. Since CRCM5-LE shows these links in 16 and 30 members, respectively, it is likely not observed in CRCM5/ERA due to potential errors in reference causal graph reconstruction (Galytska et al 2022) or internal variability. Exact agreement between reference data and model network cannot be expected due to internal climate variability. However, larger ensemble sizes increase the likelihood of capturing observed relationships within the ensemble range (compare Karmouche et al 2022).

Due to algorithm construction, the absence of links is a more robust result than their existence: if no statistical link can be established among two nodes, it is unlikely that a physical process links both under the assumption of data faithfully representing the underlying physical processes (Runge et al 2019a). All links are highly invariant towards leaving out arbitrary nodes in calculation (except for nodes contained in a given link, figure 3(b)).

We next analyze pairs of core regions before and after transitions, i.e. the links derived in figure 3. All selected transitions occur most often during July and August (figure 4(a)); however, transitions in Western and North-Western regions tend to occur earlier (beginning of July) than South-Eastern regions (beginning of August) with the exception of transition IP–WE2. This is connected to the seasonal preference of heatwave clusters themselves (Felsche et al 2022).

Zoom In Zoom Out Reset image size

The links found in figure 3 are next evaluated against all possible links between the first and second cluster in a transition (k₁ and k₂). All of them rank among the top three most frequent transitions from their respective starting cluster (figures 4(b) and (c)). Transition occurrences of IP–WE2, WE2–BI and SEE–GSI, i.e. comparably short distances, exceed the expected value of randomly occurring transitions (i.e. k₁ and k₂ being independent) by factor 3 or 4. Among these, SEE–GSI mirrors source region–target region directions of trajectories found in Bieli et al 2015. The authors also identify trajectory source regions for heatwaves in the BI region west to Great Britain, which is not represented in figure 3. Generally, source regions (i.e. start regions in our analysis) of air masses are located close to heatwave occurrence regions (Bieli et al 2015).

Expressing transition strength in terms of conditional frequencies, figure 4(c) supports the results from the causal discovery network: Values of 1 would indicate a strict connection between both clusters, i.e. k₁ is always followed by k₂ and k₂ being always preceded by k₁. In other words, heatwave occurrence in k₁ is likely a sufficient and necessary condition for heatwave occurrence in k₂. Figure 4(c) shows high cross-ensemble values (and thus suggests strong connections) for, e.g. transitions SEE–GSI or IP–WE2, but also large variability among members: In 45% of all SEE occurrences, the heatwave propagates to GSI. Equally close connections characterize transitions WE2–BI (43%) or GSI–SEE (42%). This indicates a preference of heatwaves moving from the given first core region to the corresponding second core region within the analysis period.

We next seek to connect these previously discussed transitions with potential drivers explaining their pathways (figure 5). Since negative soil moisture anomalies prevail before and after the transitions in the core regions (figure S3), we focus on centroid positions of positive Z500 anomalies during the two days preceding and two days following an indicated transition (figure 5). Our hypothesis states that a moving center of high pressure anomalies, as strong predictors of heat (Suarez-Gutierrez et al 2020), draws heatwaves with it in a similar direction. Acknowledging the fact that single events may deviate from these composites and the link between composites and heatwave occurrence may be not specific (Boschat et al 2016, Clemesha et al 2018, Zschenderlein et al 2019), we only examine regions where the signal exceeds the standard error among composite components. As was found for stationary heatwave patterns (Stefanon et al 2012), Z500 anomaly composites and temperature patterns correspond well in their spatial extent during transitions. Some heatwave footprints, however, are observed to be larger than their associated synoptic patterns (Spensberger et al 2020). The IP–WE2 Z500 and heat anomaly patterns as well as the propagation tracks are similar to the 'European Cluster' in Sánchez-Benítez et al 2020.

Zoom In Zoom Out Reset image size

In general, Z500 anomalies are located to the North-East of the heatwave core regions (compare crosses and triangles in figure 5). Blocking is known to intensify hot extremes southwest of its occurrence by various heating processes and advecting warm continental air masses (Pfahl 2014). For BI–CEE, CEE–SEE, BI–WE1, WE1–SEE, and SEE–GSI, centroid positions of Z500 anomalies lie between start and end positions of the transitions, in parts clearly mirroring the movement direction (see positions of 1–4 relative to each other, figure 5, and figure S3 for a comparison of selected links with CRCM5/ERA).

These associated patterns tend to appear during 12%–33% of the respective transitions in CRCM5-LE. Unsurprisingly, relative occurrences of the indicated transitions during all pattern occurrences remain low (0.1%–0.8%), given a high frequency of these patterns during May–October. Thus, the specific Z500 pattern occurrence is very likely not a sufficient cause for transition occurrence, but may be treated as a necessary condition (see Boschat et al 2016, Hannart et al 2016). Among pattern occurrence during any heatwave transition, however, the indicated transitions occur up to ten times more frequently, ranging from 2% to 7%. To summarize, associated pattern occurrence alone cannot be interpreted as a clear sufficient or necessary condition to causally explain heatwave transitions, although they tend to propagate in the same direction as contemporaneous heatwaves.