Robust intra-model teleconnection patterns for extreme heatwaves

Introduction: The goal of this study is to provide analysis of statistics and dynamics of extreme heatwaves over two areas of Europe, France and Scandinavia, while comparing and contrasting the representation in climate models and reanalysis.

Methods: The 1000 year long datasets are generated using respectively two climate models of different complexity. The composite maps of the heatwaves found in these datasets are compared to the ones observed in reanalysis by estimating significance of such patterns. We also employ time series analysis, in particular compare return time plots and use Gaussian stochastic processes to model the temporal correlations for rare events.

Results: Our investigation reveals that recurrent wavenumber three teleconnection patterns distinctly underpin heatwaves in France and Scandinavia. These patterns manifest in both surface temperature and mid-tropospheric geopotential height. For heatwaves of return times of 4 years such patterns are robustly represented across diverse models of complexity and reanalysis data. For longer return times, reanalysis records are too short to give statistically significant results, while models confirm the relevance of these large scale patterns for the most extreme heatwaves. Moreover, A time series analysis shows that heatwave indices defined at synoptic scale are fairly well described by Gaussian stochastic processes, and that these Gaussian processes reproduce well return time plots even for very rare events.

Discussion: These results suggest that extreme heatwaves over different areas of Europe show recurrent typical behaviours in terms of long-range spatial correlations and subseasonal-scale temporal correlations. These properties are consistently represented among models of different complexity and observations, thus suggesting their relevance for a better understanding of the drivers and causes of the occurrence of extreme midlatitude heatwaves and their predictability.

1 Introduction

Some of the most severe impacts of climate change are caused by rare and extreme events. Increased frequency and magnitude of extreme heatwaves is one of the most immediate and significant effects of global warming Van Oldenborgh et al. (2022). In the last decades a number of record breaking heatwaves have been observed IPCC (2014); Seneviratne et al. (2021). In the Northern Hemisphere mid-latitudes examples include the Western European heatwave of summer 2003, with a death toll of about 70, 000 García-Herrera et al. (2010), the mega-heatwave over Russia of summer 2010 Otto et al. (2012) with the deathtoll of about 100,000, and more recently the heatwave 2021 Western North America heatwave Philip et al. (2021).

Recent studies have identified increased frequency of heatwaves occurring in mid-latitudes in Europe of duration longer than 6 consecutive days Rousi et al. (2022). In fact, the reasons for the high impact of the 2003 Western European heatwave was not only the magnitude of the temperature anomalies, but also its long duration (two successive heat events along an overall period of 1 month). In the case of the 2010 Russian heatwave, there was the compounding effect of high temperature, long duration (1 month), and related wildfires. 2021 Western North America heatwave lasted about 10 days with peculiarly strong peak on a 2-day scale. In this case the main cause of the impacts where short term extreme temperature fluctuations, that were unprecedented in historical record World Weather Attribution (2021); Schiermeier (2021).

Heatwaves typically occur when a stationary high pressure anomaly over a region leads to subsidence and increased incoming shortwave radiation fluxes, that increase surface and near surface temperatures (Perkins, 2015; Horton et al., 2016). Short-term events lasting a few days are typically understood as caused by the occurrence of an atmospheric blocking over the region impacted by the heatwave. The most extreme and persistent events, on the other hand, require the activation of more particular dynamics and of feedback processes that enhance the surface temperature response and the duration of the event. This involves processes acting on space and time scales that go beyond the regional spatial scales and the synoptic time scales characteristic of standard heatwaves. A comprehensive analysis of the properties of extreme midlatitude heatwaves however is hindered by their rarity, the inherent noisiness of midlatitude weather fluctuations, and thus the difficulty to sample a sufficient number of events to provide a reliable statistics.

From a spatial point of view, it has been observed in the literature in the past years that some of the most extreme and persistent heatwaves seem to be associated with large scale atmospheric teleconnection patterns over the entire Northern hemisphere, for example, extreme heatwaves over the central USA Teng et al. (2013), Alberta Petoukhov et al. (2018), and Western Europe Kornhuber et al. (2019). Following a recent classification of compound extreme events Zscheischler et al. (2020), this would correspond to spatially compounding events. Some authors have interpreted the presence of these patterns in terms of amplification of quasi-stationary Rossby waves, which are claimed to be related to the genesis of extreme events in different regions of the world during the same season Lau and Kim (2012); Petoukhov et al. (2013), Petoukhov et al. (2016); Schubert et al. (2011).

The research in this area has been directed towards comprehending the pivotal role of waves that exhibit a wavenumber ranging from 5 to 8. To this end, an array of diverse detection techniques has been utilized encompassing empirical orthogonal functions (EOFs) Teng et al. (2013), spectral analysis Kornhuber et al. (2019); Petoukhov et al. (2018), and indicators that rely on resonance models Petoukhov et al. (2013), Petoukhov et al. (2016). Despite these strides, the precise identification of these teleconnection patterns is hindered by the limited amount of events available. Furthermore, we are increasingly facing the questions ranging from estimating risks for similar “black swan“ events found in the tails of the distribution: from estimating their return times to attributing them to climate change Philip et al. (2020). This highlights the need to devise more efficient tools for extreme return time estimation.

To overcome the sampling issue, Ragone et al. (2018); Ragone and Bouchet (2020), Ragone and Bouchet (2021) have adopted rare event algorithm applied to simulations with climate models of different complexity. In this way they found that the warmest summers and most persistent heatwaves are robustly associated with teleconnection patterns with wavenumber 3-4, and were able to show that these patterns are indeed statistically significant (see in particular Ragone and Bouchet, 2021). However, a detailed analysis of the dynamical properties of these peculiar atmospheric states, and the possible relation with the amplification of higher wavenumber Rossby waves suggested in the literature, have not been investigated yet.

From a temporal point of view, there are two factors that characterise the most extreme and persistent events, up to the seasonal scale. First, these events typically occur when a train of heatwaves leads to an extreme hot summer, as in western Europe in 2003 Cassou et al. (2005). This would correspond to temporal compounding mechanisms in Zscheischler et al. (2020). Possible reasons why a given season could be more likely to develop multiple heatwaves are not clear, but they would be of great interest for what concerns the seasonal to decadal predictability of temperature extremes. Second, it is known that pre-existing anomalously dry soil conditions can amplify the intensity of heatwaves occurring at a later time, a typical mechanism, for example, of certain classes of heatwaves in Europe Stéfanon et al. (2012). This would correspond to preconditioning mechanisms in Zscheischler et al. (2020). Preconditioning is typically studied with ad hoc numerical experiments initialized with dry or moist soil initial conditions, or in terms of correlations. A more detailed analysis of the properties of the decorrelation function of regional surface temperatures over different time scales seem to be absent from the literature.

The goal of this paper is to better understand the properties of long range correlations in space and time characterizing extreme heatwaves, and to compare their robustness among models of different complexity and reanalysis data. We adopt a very simple definition of a heatwave event that specifically aims at the properties of time persistence of the temperature anomalies from the subseasonal to seasonal scale (S2S). We analyse events over France and Scandinavia. In Section 2 we describe the data, the heatwave index, and the statistical tests and analysis adopted in this work. In Section 3 we describe the teleconnection patterns we find associated with the extreme events, we show the properties of the atmospheric activity associated with these events in terms of Hayashi spectra, and we discuss the statistical significance of the patterns. In Section 4 we analyse the temporal evolution of the events and the properties of the time autocovariance function of the regional temperatures, comparing them against a first order autoregressive process. Finally in Section 5 we discuss our conclusions.

2 Methodology

2.1 Data

We study extremely rare long lasting heatwaves over different areas of Europe during the June-July-August period (JJA). The data consists of daily averages of 2 m temperature and 500 hPa geopotential height, from three different datasets: one intermediate complexity climate model (PlaSim), one CMIP5-type Earth system model (CESM 1.2), and the ERA5 reanalysis.

2.1.1 PlaSim

The Planet Simulator Fraedrich et al. (2005), Fraedrich et al. (1998), also known as PlaSim, is a climate model of intermediate complexity that has been developed to simulate the Earth’s climate. The model’s dynamical core uses a spectral transform method to solve the primitive equations governing vorticity, divergence, temperature, and surface pressure. The model has a horizontal resolution of T42 in spectral space, which translates to a spatial resolution of 2.8° by 2.8° or to 64 × 128 grid-points. Additionally, the model includes 10 vertical layers and incorporates simple parameterizations of the most important physical processes that influence the climate system, such as large-scale precipitation, clouds, moist and dry convection, boundary layer fluxes of latent and sensible heat, and vertical and horizontal diffusion and radiation. The atmospheric model is coupled to a simple bucket land surface scheme. To produce a stationary state climate close to that of the 1990s, the model is driven by prescribed values of sea surface temperatures, sea ice cover, greenhouse gas concentration, and incoming solar radiation fluxes at the top of the atmosphere. We used the model to generate 1,000 years of data in this stationary state, this dataset has also been utilized in previous studies conducted by Ragone et al. (2018) and Jacques-Dumas et al. (2022).

2.1.2 Community earth system model (CESM)

We use control runs from Ragone and Bouchet (2021) which were performed with version 1.2.2 of the Community Earth System Model (CESM) developed by Hurrell et al. (2013) utilizing an atmosphere and land only setup, with the active components being version 4 of the Community Atmospheric Model (CAM4) and version 2 of the Community Land Model (CLM2). The model was run in a statistically stationary state, with the concentration of greenhouse gases, sea surface temperature (SST), and sea ice cover held at values corresponding to the climate of year 2000. The model features a horizontal resolution of 0.9 and 1.25° in latitude and longitude, respectively, and includes 26 vertical layers in hybrid pressure coordinates.

2.1.3 ERA5 reanalysis

We consider data from the ERA5 reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) Hersbach et al. (2020); Muñoz Sabater et al. (2021), from 1950 to 2020. The data was downloaded from https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5 and re-gridded to a lower resolution of 241 by 480 while performing daily averaging over available 3 h periods.

2.2 Definition of heatwaves and return times

To define heatwaves, several indices have been used in the literature, for different purposes Perkins (2015). Many meteorological criteria used in climate studies of temperature extremes focus on sub-daily fluctuations (see for instance IPCC (2014)). However, long-lasting heatwaves are the most detrimental to health and biodiversity Barriopedro et al. (2011). Moreover, many of the extreme heatwaves with the largest impact, for instance, the Western European one in 2003 or the Russian one in 2010 lasted long, from two to 5 weeks. These long lasting heatwaves were often composed of several sub-events compatible with the classical, short-time based definitions Perkins (2015). The lack of comprehensive studies of the statistics of long-lasting events has actually been stressed in the last IPCC report Seneviratne et al. (2021). Moreover, many definitions that actually involve a measure related to the persistence of anomalous daily maximum temperature values with prescribed amplitude do not always carry a natural definition of a heatwave amplitude Perkins (2015). This prevents to study independently impact of amplitude and duration of the heatwave and calls for a complementary definition of heatwaves, that can quantify both their amplitude in terms of temperature and their duration, in an independent way.

In this study we choose a specific criterion for selecting heatwaves (following the definitions of Ragone et al. (2018); Gálfi et al. (2019); Ragone and Bouchet (2020); Gálfi and Lucarini (2021); Ragone and Bouchet (2021); Gálfi et al. (2021)), which consists of extremes of time averaged surface temperature fluctuations (anomalies) defined as

$A\left(t\right)≔\frac{1}{T}{\int }_t^{t+T}\frac{1}{\left|\mathcal{D}\right|}{\int }_D\left(T_s-\mathbb{E}\left(T_s\right)\right)\left(\vec{r},u\right)\mathrm{d}\vec{r}\mathrm{d}u(1)$

where $\mathbb{E}\left(T_s\right)(\vec{r},t)$ is the mean surface temperature at each point of the grid, with a seasonal and spatial variation. In the following the symbol $\mathbb{E}$ refers to the statistical empirical average over all available years, which corresponds to the duration of the datasets (1,000 years for PlaSim and CESM and 71 years for ERA5. In this study the heatwave area $\mathcal{D}$ correspond either to France or Scandinavia and is depicted via shaded boxes on Figure 9. France is defined as a land enclosed between 43 N and 51 N and 4 W–6 E, while Scandinavia is defined as a land enclosed between 5 E − 39 E and 57 N–71 N. The heatwave duration is chosen depending on the impact of interest. In this study we select T = {5, 14, 30, 90} days. In order to remove the climate change signal from ERA5 reanalysis the fields were linearly detrended for each grid point and anomalies were defined with respect to such detrended mean.

The choice to study the extremes of temperature anomalies rather absolute temperature is made for the following reason. From a dynamical perspective, the summer months are similar, and provide more statistics than individual months. This definition emphasizes dynamical characteristics, which are thought to be described reasonably well on a seasonal time scale, rather than physical impacts which are often related to the physical temperature. Moreover, the emphasis on anomalies allows to collect richer statistics and to concentrate on properties that are better captured by the models, as opposed to the absolute magnitudes.

We are interested in summer statistics June July August (JJA) so that duration of heatwave is constrained to JJA and cannot occur outside this period. The yearly block-maxima extremes a_i occurring at calendar day t_i are defined as

$a_i≔A\left(t_i\right)=\underset{t\in \text{summer\hspace{0.17em}}\left(i\right)}{\max }\left\{A\left(t\right)\right\},(2)$

Since temperatures are correlated on a daily timescales one has to be careful about the selection of the maxima which occur on the first or the last summer day. Thus, for such events we check whether they correspond to local maxima. If they do not, we look for greatest local maximum in the summer time series for that year, thus avoiding registering extrema of May or September.¹

To compute the return times, the yearly summer extremes a_i are ranked in decreasing order, as usual (see Section 4.2.1). Next, we compute return times for heatwaves using the method Lestang et al. (2018) described in Supplementary Material. The return time expression reads in our notation as

$r=\frac{1}{\log \left(1-\text{\hspace{0.17em}rank}\left(a\right)/M\right)},(3)$

where M corresponds to the total number of years and rank(a) corresponds to the index i of the extreme found in the sequence $\left\{a_i\right\}$. For very rare values a, this expression gives r ≈ M/rank(a), which is the more familiar definition of return time found in the climate literature.

2.3 Statistical analysis

To identify long lasting heatwaves we will consider time averages of temperature and geopotential height on different time scales. In particular, given a variable X(x, t) we consider averages conditioned on the extreme, also called composite maps. The composite map of X at the lead time τ of an extreme heatwave with a return time r is

$X^r\left(\mathbf{x},\tau \right)=\frac{1}{N_r}{\displaystyle \sum _{\left\{i\left|a_i\ge a^r\right.\right\}}}X\left(\mathbf{x},t_i-\tau \right).(4)$

where N_r is a number of events for a given threshold a_r and t_i is the time of the onset of the heatwave. This definition gives informations on how a typical heatwave above a given threshold develops and disappears.

We test statistical significance using the standard Student t test. For instance, let us discuss the estimation of the average of a random variable Y, using N independent samples {Y_n} through the empirical average

${\mathbb{E}}_e\left(Y\right)=\frac{1}{N}{\displaystyle \sum _{n=1}^N}{Y}_n(5)$

If our aim is to test whether this sample average is significant, in the sense that the probability that the actual average is 0 is very low, we will compute the t-value

$t=\frac{\sqrt{N}{\mathbb{E}}_e\left(Y\right)}{{\sigma }_e\left(Y\right)},(6)$

where σ_e(Y) is the empirical standard deviation

${\sigma }_e\left(Y\right)=\sqrt{\frac{1}{N-1}{\displaystyle \sum _{n=1}^N}{\left[Y_n-{\mathbb{E}}_e\left(Y\right)\right]}^2}.(7)$

We consider the average statistically significant with |t|≥ 2. We note here that we have taken a threshold on the t-value rather than on the confidence probability. However our smallest sample size is around 16, and the t-value for a given confidence probability changes very weakly with the sample size above 15–20. In any case with our sample sizes a t-value of 2 is conservatively consistent with 95% confidence level at the very least (Hans von Storch, 2002).

The method of estimating uncertainties known as bootstrapping will be used. This involves breaking the total temporal sequence into 10 subsets and computing the relevant quantities correspondingly. Afterwards, mean ⟨x⟩ serves as the best estimate, while standard deviation σ_x as uncertainty.

3 Hemispheric teleconnection patterns

3.1 Hayashi spectra of midlatitude atmospheric waves

The extreme teleconnection patterns found by Ragone et al. (2018); Ragone and Bouchet (2020), Ragone and Bouchet (2021) in simulations with climate models of different complexity, suggest that the local atmospheric blockings responsible for heatwaves occurring in subregions of Europe could be associated with stationary or very slowly moving atmospheric waves with wavenumber between 3 and 4. The body of literature of Petoukhov et al. (2013), Petoukhov et al. (2016), Petoukhov et al. (2018); Kornhuber et al. (2019), Kornhuber et al. (2020) identifies the amplification of quasi-stationary Rossby waves with a different range of wavenumbers, between 5 and 8 as responsible for the occurrence of several observed midlatitude extreme events, including heatwaves.

Since we are seeking to detect global quasi-stationary Rossby wave patterns we use classical space-time Fourier decomposition into eastward/westward and stationary Rossby waves Hayashi (1971). This analysis is performed on 500 hPa gepotential height averaged in the meridional direction over a latitudinal band centered over the region affected by the heatwaves. The procedure is done via an assumption about the nature of the waves, such as attributing incoherent part of the spectrum to real travelling waves Pratt (1976); Hayashi (1979). The low frequency - low wavenumber domain corresponds to standing and westward propagating waves, while the high frequency-high wavenumber to synoptic disturbances. For the details also see the appendix of Dell’Aquila et al. (2005). We chose the latitudinal belt 55–75 N as a region relevant for Scandinavia over which the averaging was performed. This belt is part of the larger belt 30–75 N where most of the baroclinic and the low frequency wave activity occurs. Following the procedures in the literature, the meridional integration is an average, not a sum, so the total power does not scale with the size of the latitude band. Additionally, in order to compensate for the non-constant density of points on a log-log plot the spectra are multiplied by k ⋅ ω/2π.

In what follows we work with conditional averages (condition is to keep only the years where A(t) > th, Eq. 1, for a certain threshold th, so that th = −∞ corresponds to climatology). We show Hayashi spectra for different thresholds on Figure 1. On the left panel we plot climatology threshold th = −∞, which displays the classical spectrum of Rossby waves. On the right panel we plot Hayashi spectra for the summers with 16 most extreme T = 30 days heatwaves, which corresponds to the threshold th = 4.5 K. We observe a structure at a period of order 50 days and wavenumber 3 for eastward propagating waves. We note that this pattern in the Hayashi spectrum is isolated from the rest of the Rossby wave branch and it suggests that quasi-stationary processes take place when conditioned to long-lasting heatwaves.

FIGURE 1

FIGURE 1. Comparison of Hayashi spectra for (A) 1,000 years of CESM climatology and (B) a subset of 1,000 years of CESM with Scandinavian heatwaves above the threshold of 4.5 degrees of T = 30 days heatwaves. Boundaries of the latitude band over which the meridional average is taken are 55 and 75.

3.2 Extreme teleconnection pattern

Hayashi spectra associated with heatwaves in Scandinavia, computed with the CESM model, show a clear local, low-frequency, low wave-number peak (See Figure 1). There are two questions that arise immediately.

• What is the real-space pattern that corresponds to this peak?

• Does the teleconnection pattern have a counterpart in the observational record

The chief quantities of interest are 500 hPa geopotential height (500GPH) and the 2 m temperature (T2M). The 500GPH is a convenient quantity to analyse the large scale circulation in the troposphere where the effects of the topography are minor Blackmon (1976).

To answer the first question we compute composite maps (see Equation 4) of 500GPH and T2M on Figure 2 conditioned to the T = 30 days long heatwaves of Scandinavia of threshold a = 3.5 K (containing 65 events). The choice of the threshold is motivated by the comparisons with the actual event which was observed in reanalysis that will be commented on below. We observe a tripole-like structure, which consists of an Arctic cyclonic anomaly surrounded by three anti-cyclonic anomalies, over Scandinavia, Chukotka and North-East Canada. The cyclonic anomaly extends to Greenland and Siberia, and it is associated with negative temperature anomalies. This teleconnection pattern indeed corresponds to the local maximum in the low-frequency low-wavenumber branch of the Hayashi spectrum (Figure 1B).

FIGURE 2

FIGURE 2. (A) Northern Hemisphere surface temperature anomaly (colors) and 500-hPa geopotential height anomaly (contours), conditional on the occurrence of European heatwaves with T =90 days and the threshold of a =2 K, estimated from the large deviation algorithm, courtesy to Ragone et al. (2018). (B) ERA5 reanalysis averaged over July 2018 and detrended relative to the 1950–2020 linear trend, we follow the same conventions for representing the geopotential height and temperature anomalies, except that negative geopotential anomalies are presented via dashed lines. (C) Composite map of CESM drawn from 1,000 years long run and conditioned to T =30 days heatwaves above the threshold of 3.5°. We follow the same conventions for representing the geopotential height and temperature anomalies. The temperature colormap is shared among the three plots.

Motivated by the second question we select the July 2018 event from ERA5 reanalysis. The details concerning the data processing of ERA5 are discussed in the Section 2.1.3. The July 2018 event consists of a prolonged heatwave in Scandinavia, whose temperature anomaly was above 3–3.5° for about a month. Its composite map (Figure 2B) displays various similarities with the CESM teleconnection (Figure 2C). Most striking resemblance are the negative Arctic 500GPH anomaly, as well as positive anomalies in Scandinavia and North-East Canada. Reanalysis favors a positive 500GPHA shifted eastwards towards the Atlantic. In the Pacific ocean we see a train of 500GPHA instead of a single positive anomaly. In general, one finds that geopotential anomalies are stronger in reanalysis; for instance, 500GPHA over Scandinavia is about 140m, whereas in CESM we find values on the order of 70m. The stronger anomalies can be attributed to the fact that in reanalysis we were averaging over the single event, while in CESM the averaging is performed over multiple heatwaves.

We are now in a position to also comment on the relative similarity between these two teleconnection patterns with the one presented in Ragone et al. (2018). For illustration purposes we display the copy of the relevant teleconnection on Figure 2A. It was obtained using the rare event algorithm for PlaSim conditioned to T = 90 days heatwaves of threshold a = 2 K occurring in Europe (to be more precise each grid point is averaged over 90 days). It turns out that the main contribution is due to heatwaves in Scandinavia, of threshold that can be estimated as a ∼ 3.5–4 K. It is interesting that the similarities between the panels exist despite the different complexity of the two models, PlaSim and CESM. This gives us indication that such teleconnection patterns could be robust feature of the large scale dynamics. However, the specific conditions are not identical, and the similarity to the July 2018 event could be accidental. For instance, the high temperature and anti-cyclonic anomaly in the southeastern Asia is consistently reproduced in PlaSim and ERA5 but not in CESM. In order to better test the assertion that planetary-scale teleconnections are consistent, we perform a more extensive intra-model comparison in the next Section.

3.3 Teleconnections in models versus reanalysis

We aim to provide systematic comparisons between the models. We concentrate on two areas in Europe, corresponding to France and Scandinavia. The areas have been chosen to correspond closely to the heatwave clusters in Western Europe and Scandinavia identified by Stéfanon et al. (2012).

Since the ERA5 dataset contains only 71 years we limit the study to return times not larger than 4 years, in order to have enough data (at least 18 heatwaves). For consistency we select heatwaves with the same return time of 4 years in CESM and PlaSim (250 heatwave events in both models) which allows for stronger statistical inference. The corresponding thresholds a are not the same across the models, especially in Scandinavia, where the threshold for PlaSim is much lower than the other two datasets (see Table 1). This is further discussed in Section 4.2.1.

TABLE 1

TABLE 1. Values of a r =4 years return time threshold a₄ computed for France and Scandinavian heatwaves (columns) over PlaSim, CESM and ERA5 datasets (raws).

One wonders how similar are the teleconnection patterns drawn from models at the onset of individual heatwave events, i.e., at τ = 0 days. We take duration T = 14 days, and plot the corresponding composites on Figure 3, following the method described in Section 2, for heatwaves occurring in France (top) and Scandinavia (bottom). We clearly see the presence of a tripole-like structure in all the images, but there are many additional details. PlaSim and CESM share more features, although qualitative similarities can also been seen when comparing with ERA5. In particular, one observes statistically significant positive 500GPHA in the North sea of order 84 m in PlaSim, 90 m in CESM and 50 m in ERA5. The second largest in magnitude feature is the negative Arctic 500GPHA of order −50 m in PlaSim and −82 m in CESM which extends to the Atlantic. The extension of this pattern seems to be consistent with what one observes in ERA5. However, in ERA5 it is not statistically significant. In contrast, we have a statistically significant negative 500GPHA in the Atlantic ocean of −133 m, which is not so strong in PlaSim or CESM. The third order feature is the positive 500GPHA in East Canada which reaches 30 m in CESM and PlaSim and approximately 40 m in ERA5. Finally, in PlaSim we have 10 m in North-East Asia anomalies, which are stronger in CESM and ERA5. Most features are statistically significant for CESM and PlaSim but many are not for ERA5 (see Section 2.3 on statistical significance test), which is due to scarcity of the observational record.

FIGURE 3

FIGURE 3. Composites of r =4 years return time in different models (left to right), where 2 m temperature and 500 mbar geopotential anomalies are plotted conditioned to different areas (top-bottom): (A) PlaSim France. (B) CESM France. (C) ERA5 France. (D) PlaSim Scandinavia. (E) CESM Scandinavia. (F) ERA5 Scandinavia. at τ =0 lag (onset of the heatwave). The colormap for the temperature is indicated next to the figures. Only |t|>2 statistically significant anomalies of temperature are displayed (see the colormap in K), whereas all level sets of geopotential height anomalies graded by 10m are shown (Green - positive, Orange - Negative). The values of geopotential that are statistically significant |t| > 2 are plotted with solid lines and otherwise with dashed lines (see Section 2.3 on statistical significance). The global maximum and minimum of geopotential is shown using blue and red text tooltips respectively. All subsequent figures displaying composites will follow this general style.

The panels at the bottom of Figure 3 corresponding to teleconneciton patterns conditioned to Scandinavian heatwaves also show of these three features: a primary one consisting of positive 500GPHA in Scandinavia: 96 m in PlaSim, 110 m in Scandinavia and 114 m in ERA5; a secondary one consisting of negative 500GPHA in the Arctic with −39 m in PlaSim, −73 m in CESM, and in ERA5 of order—90 m. While this peak is larger than expected, we find that overall shape and magnitude of 500GPHA in the Arctic are qualitatively consistent between the models and reanalysis. The tertiary pattern consists of positive 500GPHA in Eastern Canada which is approximately 20 m in both PlaSim and CESM and of order 60 m in ERA5 with a easterly shift. Finally there is also a pattern which seems to be generally consistent across the models and reanalysis: the Mediterranean depression of order −20 m. However, in reanalysis it is, again, not statistically significant.

To summarize, we find qualitative agreement between Figures 3A–C as well as agreement between 3d, 3e, and 3f. CESM seems to be capturing the teleconnection patterns closer to reanalysis. For instance, in case of heatwaves conditioned to France we see that CESM captures better the distribution of geopotential anomalies over Canada and Northern Russia when compared to ERA5. Nevertheless, the similarities displayed in the teleconnection patterns suggest that they are model independent and have observational counterparts that are rather robust.

Another important question is how these patterns change when conditioning for more extreme heatwaves. In Figure 4 we compare extreme heatwaves with return time larger than 100 years, with the ones with return time in the range 4 < r < 100. In the latter case we have actually removed the 10 most extreme events to show that they do not contribute noticeably compared to the remaining 240 events. Like in the previous analysis we identify the three most important features as negative GPHA over Greenland and positive GPHA over North sea and East Canada. In Figure 4A, for the events that are more extreme, they take values 166 m, −134 m and approximately 90 m, while on Figure 4B, for the events that are milder, they take the values 87 m, −81 m and approximately 40 m. Furthermore, Figure 4A has variety of other mostly statistically insignificant cyclonic/anticyclonic anomalies. The conclusion is that when going from less extreme to more extreme heatwaves, the most significant features remain geographically fixed (500GPHA naturally deepen), while we obtain a few more anomalies that do not pass |t| > 2 test. For comparisons with Scandinavia see Supplementary Figure S1.

FIGURE 4

FIGURE 4. Comparison between (A) r >100 and (B) 4< r <100 for France heatwaves in CESM dataset. See the caption of Figure 3 for plotting conventions.

Finally, we assess how the teleconnection patterns depend on the choice of the heatwave duration T (see Eq. 1). We test 4 periods, T = {1, 5, 14, 30} days for r = 10 years return times, that are shown in Figure 5 for the CESM dataset conditioned to heatwaves in France. The patterns discussed before, including the nodes of the tripole structure, such as positive GPHA on the North sea, Canada and East Siberia can be recovered also from Figures 5B,D with mostly the same intensities and some minor shifts. Most different of all is Figure 5A which is actually a daily average. In this case the Arctic anomaly descends towards Scandinavia/western Russia region and we have an isolated negative GPHA in the Atlantic. The differences are to be expected since on daily time scale we may observe many of the synoptic features. We conclude that the features we have described in CESM at τ = 0 appear to be also robust across different time-averaging intervals, namely, at T = {5, 14, 30} days. However, we observe that for shorter periods such as 1–5 days negative anomaly in the Arctic develops tongues that reach the midlatitudes. These are not present for large values of T as the faster synoptic type perturbations are averaged over. For comparisons with Scandinavia see Supplementary Figure S2.

FIGURE 5

FIGURE 5. Composites of 2 m temperature and 500 mbar geopotential conditioned on r ≥10 years return time heatwaves over France at τ =0 lag. The reference model is CESM. Among the different panels we vary the length of heatwaves: (A) T =1, (B) T =5, (C) T =14 and (D) T =30 days. See the caption of Figure 3 for further plotting conventions.

4 Temporal evolution of the extremes

The goal of this Section is to compare the temporal statistics between the reanalysis and the models. The time evolution of composite statistics will only be considered for CESM, whereas other types of statistical properties will be compared across the models.

4.1 The dynamical evolution of composites

Here we address the dynamical evolution leading up to and during the heatwave. We limit this study to the CESM model because we lack data in ERA5 reanalysis and CESM is a higher fidelity model than PlaSim. In Figure 6 we show the dynamical evolution of heatwave composites, conditioned to r ≥ 4, in France, and varying the lag time τ. This shows how the teleconnection patterns change as a function of lead time. We see that 15 days prior to the onset of heatwaves (Figure 6A) we already have significant negative GPHA of order −27 m in the Arctic region that slowly broadens and deepens (−35m at τ = −6 days, Figure 6B, and −67 m at τ = −3 days; Figure 6C). Meanwhile, an anticyclonic block develops at the midlatitudes, reminiscent of a similar pattern in Figure 5C. Five days after the onset (Figure 6D) we can distinguish the tripole-structure (see Section 3.2) with strong significant 133 m positive GPHA in the North sea and negative GPHA that attains −76 m in the North Pole. This pattern starts to disappear at τ = 15 days (Figure 6F), corresponding to the duration of the heatwaves.

FIGURE 6

FIGURE 6. Composites of 2 m temperature and 500 mbar geopotential conditioned on r ≥4 years return time heatwaves in France at τ =0 lag. The reference model is CESM and T =14 days is chosen. Control parameter for this figure is τ with (A) τ =−15, (B) τ =−6, (C) τ =−3, (D) τ =5, (E) τ =10 and (F) τ =15 days, thus time evolution of CESM composite heatwave in France is displayed. See the caption of Figure 3 for plotting conventions.

4.2 Time series analysis

We discuss the properties of the time series of the area averaged temperature, over France and Scandinavia separately. The idea is to show how closely the models can resemble the real time series and, further, whether we can reduce the dynamics to a simple stochastic process.

4.2.1 Return time plots

We show the probability distributions of the 14 days running mean of 2 m temperature anomaly (T2MA) over France and Scandinavia on Figure 7. Both areas have Probability Distribution Functions (PDFs) that are Gaussian in the bulk, with fatter tails, particularly for PlaSim (see m Section 4.2.3 for more details on this). In general, reanalysis datasets are characterized by a noticeable skewness but relatively small kurtosis, which is captured well by CESM in both areas. PlaSim, on the other hand, does not represent well PDF of ERA5 in Scandinavia. Nevertheless, Gaussian statistics seems a reasonable initial assumption as we shall see in Section 4.2.3. This is consistent with the fact that we perform long-time large-area averages: time-averaged statistics that are long enough should assymptote to Gaussian distribution, assuming they are not too correlated.

FIGURE 7

FIGURE 7. 14 days running mean distributions of daily 2 m temperature anomaly (T2MA) integrated over the area of (A) France, (B) Scandinavia, in 3 different datasets: (blue) JJA 1000 years dataset of PlaSim, (green) JJA 1000 years dataset of CESM, (red) JJA 1950–2020 ERA5 reanalysis dataset. We see that France data follows roughly Gaussian statistics, while in Scandinavia PlaSim and to some extent CESM favor non-Gaussianity. The values of standard deviation (std), skewness (skew) and kurtosis (kurt) are given in the inset for each color-coded distribution.

Figure 8 depicts return time plots computed using Eq. 3. The left panels correspond to return times for France heatwaves. In particular, Figure 8A shows a remarkable agreement between PlaSim and reanalysis, despite PlaSim being a climate model of intermediate complexity. Long-lasting 90 days heatwaves are fitted even better by CESM (Figure 8C). If we look at the most extreme 14 days events (orange dots) in ERA5 reanalysis there are two apparent outliers (the most extreme being the European heatwave 2003) compared to the climate models. However, since for these two events we are in the rightmost tail of the distribution where we don’t have enough data one could argue that the climate model fits (we described) cannot be rejected solely based on these two events.

FIGURE 8

FIGURE 8. Return time plots for heatwaves, summer extremes of 2 m temperature anomaly (T2MA), in different regions (left-right) and different models (bottom-top). (A) PlaSim France, (B) PlaSim Scandinavia, (C) CESM France and (D) CESM Scandinavia. The dots correspond to ERA5 data, while the shaded regions to bootstrapping on climate model data, i.e., the data is split into 10 subsets and return times are computed using the approach described in Section 2.2. Then, mean and standard deviation are computed so that the shaded region corresponds to mean plus or minus one standard deviation. The colors and the corresponding duration of heatwaves T are indicated on the legend.

For Scandinavia the models do not compare so well with ERA5 reanalysis. For instance, Figure 8B demonstrates that PlaSim systematically underestimates the intensity of the extreme events by a large margin and across different values of T, i,e. not only the few extreme outliers. As a result, the 4 years return time has a threshold a₄ = 1.74 K in PlaSim and a₄ = 3.25 K in CESM consistent with the Table 1. Figure 8D shows that CESM performs much better for Scandinavia. The return times of the events are overestimated, at least for T < 30 days. In other words, relationship between thresholds and return times is not the same. Nevertheless, it is much closer to reanalysis, for instance for 4 years return time events the predicted threshold is a₄ = 3.79 K (Table 1).

One may wonder why models have difficulty accurately describing Scandinavia heatwaves. From the theory of single time scale Gaussian stochastic fluctuations, return times depend on the variance and decorrelation time scale of the time series Lestang et al. (2018). The picture is more complicated when dealing with time averages of signals with multiple time scales. However, generally speaking, return times are expected to scale with the variance of the underlying probability distribution (See Section 4.2.2). Figure 9 shows standard deviation of the grid-point 2 m temperature anomaly (T2MA) PDF taken over JJA period in PlaSim CESM and ERA5. It demonstrates that the models have trouble capturing precisely the variance over the region of north Russia/Scandinavia, where PlaSim has a strong negative bias while CESM a small positive bias. This goes hand-in-hand with how ERA5 return time plots are underestimated in Figure 8B and overestimated in Figure 8D. In general, according to Figure 9, PlaSim captures poorly standard deviation in higher latitudes. If we consider other parts of the European continent differences are small and, broadly speaking, CESM compares more favorably to ERA5.

FIGURE 9

FIGURE 9. T2M standard deviation over JJA in different data sets: (A) PlaSim, (B) CESM, (C) Detrended ERA5 reanalysis. Relevant areas (France and Scandinavia) that are used below for regional analysis are shaded.

4.2.2 Autocovariance function

To have a more complete understanding of the statistics of the tails of the distribution, the standard deviation is not sufficient and the decorrelation properties of the signal must be investigated. We compute the autocovariance function of the T2M time series over France and Scandinavia in PlaSim and CESM and compare them to the corresponding time series in ERA5. The autocovariance function is given by

$C\left(t\right)≔\mathbb{E}\left[A\left(X\left(t\right)\right)A\left(X\left(0\right)\right)\right],(8)$

When the process can be described by a single time scale one expects exponential scaling and the corresponding time scale can be extracted simply as

${\tau }_c≔{\int }_0^{+\infty }\frac{C\left(t\right)}{\mathbb{E}\left[A^2\left(X\left(0\right)\right)\right]}\mathrm{d}t(9)$

A similar analysis has been performed in Ragone and Bouchet (2020) on 1,000 years of data generated with PlaSim in perpetual summer conditions. In that case it was found that time series of area averaged surface temperatures show two time-scales, the fast one of order of 4 days related to synoptic variations and the slow one of the order of 30 days, which the authors relate to heat capacity of the soil and therefore the soil moisture content.

As stated above, the properties of the autocovariance function provide a contribution to the return time plot. As is known from the theory of stochastic processes Gardiner (1985) a Gaussian process can be generated knowing the parameters of the autocovariance function. We compute it using the PlaSim and CESM time series of T2M integrated over the area of France and Scandinavia and plot the results on Figure 10. The plots show that the processes involved indeed are well represented by the two time scales, which means they can be fitted by two exponential functions. The fit parameters do not precisely match the ones provided by Ragone and Bouchet (2020) which is not surprising since in that study a heatwave over much larger European region was considered in the perpetual summer regime. Another important parameter to be inferred from Figure 10 is the cross-over scale, which corresponds to the time at which the slow time scale starts to dominate the decay.

FIGURE 10

FIGURE 10. Autocovariance function across two regions (left-right) and two models (top-bottom): (A) PlaSim France, (B) PlaSim Scandinavia, (C) CESM France and (D) CESM Scandinavia. The blue dotted curves correspond to the 1,000 years-long time series, while the blue dashed lines represent the fit given by Eq. 2. To provide an estimate of uncertainty we have performed bootstrapping (splitting the time series into 10 trajectories) indicated by the yellow region which represents mean plus or minus one standard deviation. The ERA5 autocovariance function over 71 years is displayed via red dashed line and the underlying time series by red dotted line.

The presence of two time scales motivates performing exponential fit, using the following nomenclature for its parameters

$f=A_1\exp \left(-\frac{t}{{\tau }_1}\right)+A_2\exp \left(-\frac{t}{{\tau }_2}\right),(10)$

with τ₁ < τ₂ the fast and slow time-scales respectively. The importance of the slow time scale depends on the ratio A₁/A₂. The fit is applied to PlaSim/CESM/ERA5 for both France and Scandinavia T2M series and the parameters of the fit are displayed in Table 2. According to this table both PlaSim and CESM approximate the first time scale (slope). Indeed, this can be confirmed across all the panels in Figure 10A. However, Figure 10B shows that in Scandinavia PlaSim does not get the variance of the time series right since the blue and the red curve are not matching (also τ₁ is off by approximately factor 2). Another conspicuous feature is the mismatch between the plateau in Figure 10A between PlaSim and ERA5. These discrepancies translate directly into mismatch between return time plots of PlaSim and ERA5 as displayed in Figure 8B al ready discussed earlier.

TABLE 2

TABLE 2. Exponential parameters for equation (10) fit. The coefficients are given inside the table for each model (rows) and both areas: France on the left and Scandinavia on the right.

CESM fits ERA5 quite well in both France and Scandinavia (see Figures 10C,D), although in the latter case there is an over-estimation of variance discussed earlier in Section 4.2.1. Since ERA5 time series is relatively short we cannot make robust comparisons for the second time scale, which tends to be of the order of 30–80 days and, as already mentioned, is likely associated with soil-atmosphere interactions. This quasi-plateau identified in PlaSim on Figure 10A can be also identified in CESM and ERA5 on Figure 10C, however to a lesser degree. This means that France heatwaves in PlaSim have a bias towards stronger temporal persistence than in observations. CESM fits return time for France heatwaves somewhat better than PlaSim (Figures 8A,C). Thus one may wonder to what extend fitting the second time scale is indispensable for providing a better fit. What if A₂ = 0? We will address this question in the next section.

4.2.3 Gaussian process

In this Section we compare the statistics of the time series presented above to the synthetic ones generated from Ornstein-Uhlenbeck (OU) process

$\mathrm{d}{X}_t^{\alpha ,ϵ}=-\alpha X_t\mathrm{d}t+\sqrt{2ϵ}\mathrm{d}{W}_t(11)$

which can be interpreted as reformulation of a Langevin equation with correlation time τ_c = α⁻¹ Gardiner (1985). From the corresponding Fokker-Planck equation one finds a stationary probability density P_s(x) that is a Gaussian distribution with standard deviation $\sigma =\sqrt{ϵ/\alpha }$.

$P_s\left(x\right)=\sqrt{\frac{\alpha }{\pi ϵ}}\exp \left(-\frac{\alpha x^2}{ϵ}\right)(12)$

${⟨X\left(t\right)X\left(s\right)⟩}_s=\frac{ϵ}{2\alpha }\exp \left(-\alpha |t-s|\right)(13)$

The first passage time conditioned on the stationary measure Lestang et al. (2018) can be computed analytically for this model.

However, it is not so simple to compute analytically the autocovariance of the time averaged series. We will therefore adopt an empirical approach following Herbert and Bouchet (2017). We simulate the OU process based on the parameters obtained from the exponential fit on Figure 10 and prepare our Gaussian process by a direct sum of the two time series obtained from the OU integrator.

$X_t^{\alpha ,ϵ}=X_t^{{\alpha }_1,{ϵ}_1}+X_t^{{\alpha }_2,{ϵ}_2}(14)$

The parameters are chosen as α_i = 1/τ_i and $\sqrt{2{ϵ}_i}=A_i{\alpha }_i$ (see Table 2). The synthetic time series of length 10⁵ (“days”) are generated which we then separate into 10 instances of 100 years-long sequences, each of the length 100 “days” in order to mimic the summer sequences from the models. The mean and standard deviation is computed based on these 10 instances. The resulting return times (Eq. 3) are plotted on Figure 11 against the ones from General Circulation Models (GCMs) that were shown on Figure 8. We observe that the match between the Gaussian process and GCMs is relatively good, especially in case of PlaSim over France (Figure 11B): the model points are within the error bars of the bootstrapped OU simulation. This is consistent with the fact that the corresponding fluctuations in the upper tail of the distribution are rather Gaussian (see Figure 7). Both short and long duration heatwaves are captured relatively well (within the error bars), although in some cases the Gaussian process tends to produce less extreme events.

FIGURE 11

FIGURE 11. Return time plots for heatwaves in both regions for T ={5,14,30,30} days. (A) PlaSim France single time scale, (B) PlaSim France two time scales, (C) PlaSim Scandinavia two time scales, (D) CESM France single time scale, (E) CESM France two time scales and (F) CESM Scandinavia two time scales. The dots correspond to the synthetic time series obtained from OU process (11). (left) single time scale OU process with (top) Fit parameters: A =8.101, τ =18.6 days, where the former corresponds to the total variance of the PlaSim data and the latter to the corresponding integral time scale and using the same approach (bottom) Fit parameters: A =7.168, τ =7 days (center, right) two time scale OU process obtained by a direct sum, Eq. 14, with the parameters taken from the Figure 10. This figure illustrates that we need two time scales to reproduce well the return time plots. The approach works slightly worse in PlaSim Scandinavia, where likely this is due to the strong kurtosis in the model (See Figure 7).

Now we are in a position to asses the importance of the second time scale. This is particularly interesting over France where the corresponding coefficient is larger (see Table 2). To answer this question we only retained the first term in Eq. 10 and chose the coefficient A₁ → A₁ + A₂ and A₂ → 0 to match the true variance and the integral time scale for τ₁ → τ_c, τ₂ → 0, Eq. 9. The mismatch is quite pronounced for PlaSim (Figure 11A). In the case of CESM (Figure 11D), on the other hand, the differences are modest but are still noticeable as the extremes start to diverge outside of the error bars. Thus we conclude that the second time scale is less relevant in CESM yet it helps to more accurately capture the return time plot.

5 Conclusion

In this paper we presented a spatio-temporal analysis of persistent heatwaves over France and Scandinavia. One of the main findings of this study is that 2-week long summer heatwaves are associated with robust teleconnection patterns that involve a strong positive 500 mbar geopotential anomaly (500GPHA) in the region of interest, negative 500GPHA near Greenland and positive 500GPHA in North-East Canada. The patterns are qualitatively consistent across the ERA5 reanalysis and two models of different level of complexity, PlaSim and CESM.

This quasi-stationary structure appears in Hayashi spectra as a wave-number 3 peak isolated from the general eastward-propagating Rossby wave spectrum. We note that the same wave-number 3 pattern is identified when a neural network is trained to predict heatwaves using a longer dataset generated with PlaSim Miloshevich et al. (2023) (but with diurnal cycle). Similar patterns and spectral features can be found during the July 2018 heatwave over Scandinavia. In the literature there has been a discussion Petoukhov et al. (2013), Petoukhov et al. (2016), Petoukhov et al. (2018); Kornhuber et al. (2019), Kornhuber et al. (2020) around the amplification of quasi-stationary Rossby waves with a different range of wavenumbers, between 5 and 8.

The teleconnection patterns are robust to changes in a period (length) of heatwave as long as it is sufficiently large (larger than several days). The main features (wave-number 3 pattern) are consistent across less extreme (with few year return times) and more extreme (from 4 to 100 years return times).

We have also shown that the statistical properties of area averaged surface temperature in the two considered regions can be modelled using a two time-scale Ornstein-Ulhenbeck process. The length of the model runs (1,000 years) allows us to draw conclusions on the quality of the fits (estimating uncertainty) of extreme heatwaves with returns as large as 100 years. This allows us to respond positively to the question on whether correctly capturing variance and autocovariance function of the process allows to properly infer the return time curves. The return times of climate models are compared with reanalysis which reveals quite good agreement for heatwaves in France. For the heatwaves in Scandinavia return times in CESM and ERA5 diverge for periods of order 5–14 days but are generally well captured for longer events.

These results demonstrate two things. First, for the regions involved in this study, numerical climate models even of moderate complexity can reproduce the main dynamical and statistical features of persistent heatwaves at global scale. In other words, long-lasting, time averaged phenomena are not only more important from the point of view of impacts, but also to a good extent easier to study with numerical models than shorter phenomena. Second, the presence of recurrent large scale teleconnection patterns during long lasting heatwaves over these regions indicates that these events are characterized by a high degree of typicality. This could be potentially very useful in order to identify precursors for the risk of extreme persistent heatwaves, thus improving their predictability, if the dynamical reasons for the formation of these patterns could be properly identified. Recent studies have suggested that typicality of heatwave teleconnection patterns found in reanalysis and CMIP6 models may be compatible with the concept of instanton in large deviation theory Gálfi et al. (2021); Gálfi and Lucarini (2021). While this is a suggestive proposition, we stress that comprehensive quantitative analysis, beyond visual similarity of two-dimensional maps of selected fields, is necessary to support this type of hypothesis.

In this study we have compared detrended ERA5 reanalysis time series with time series from simulations performed at stationary state, with fixed CO² concentration and prescribed sea surface temperatures. It would be valuable to examine how teleconnection patterns may change under different greenhouse gas emission scenarios and how these changes could impact the frequency, intensity, and duration of heatwaves. In general, the northward shift in jet streams Osman et al. (2021) is expected to contribute to drying regimes in Europe. A substantial increase in high pressure systems over UK in late summers of the second half of the 21st century is also projected Rousi et al. (2021). This regime is linked with particularly dry conditions in western Europe and may be caused by slow-down of Atlantic Meridional Overturning Circulation (AMOC) according to some studies Haarsma et al. (2015); Duchez et al. (2016). The regime may also originate as a result of circumglobal Rossby wave patterns associated with meandering jetstream Kornhuber et al. (2019), Kornhuber et al. (2020). Finally, changes in Eurasian snow cover fraction and shrinking sea ice are also significant factors increasing the likelihood of the more persistent European blocking affects and thus affecting European heatwaves in future warming scenarios Zhang et al. (2020).

Understanding the effects of future warming on the appearance and behavior of quasi-stationary Rossby waves that induce heatwaves could shed light on the potential risks associated with such events. In order to further investigate the physical mechanisms underlying the teleconnection patterns observed in this study, an analysis of the output of CMIP6 models could be used to pinpoint the interactions between these patterns and other modes of climate variability, such as North Atlantic Oscillation (NAO), considering that coupled models seem to improve the representation of teleconnection patterns Rousi et al. (2021).

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

FB and FR contributed to conception and design of the study. FR performed numerical simulations. GM organised the database. GM and PR-P performed the statistical analysis. GM wrote the first draft of the manuscript and created the notebooks on GitHub. GM, FR, and FB wrote sections of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the ANR Grant SAMPRACE, project ANR-20-CE01-0008-01 (F. Bouchet). This work has received funding through the ACADEMICS Grant of the IDEXLYON, project of the Université de Lyon, PIA operated by ANR-16-IDEX-0005. This work was granted access to the HPC resources of CINES under the DARI allocations A0050110575, A0070110575, A0090110575, and A0110110575 made by GENCI.

Acknowledgments

We acknowledge CBP IT test platform (ENS de Lyon, France) for ML facilities and GPU devices, operating the SIDUS solution E. Quemener (2014). We acknowledge the help of Alessandro Lovo in maintaining the GitHub page.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1235579/full#supplementary-material

Footnotes

¹This is because with the original definition any extreme which actually occurs in May or September can superficially inflate the summer extremes for the first day or the last day of the summer because of temporal correlations.

References

Barriopedro, D., Fischer, E. M., Luterbacher, J., Trigo, R. M., and García-Herrera, R. (2011). The hot summer of 2010: redrawing the temperature record map of europe. Science 332, 220–224. doi:10.1126/science.1201224

PubMed Abstract | CrossRef Full Text | Google Scholar

Blackmon, M. L. (1976). A climatological spectral study of the 500 mb geopotential height of the northern hemisphere. J. Atmos. Sci. 33, 1607–1623. doi:10.1175/1520-0469(1976)033⟨1607:ACSSOT⟩2.0.CO;2

CrossRef Full Text | Google Scholar

Cassou, C., Terray, L., and Phillips, A. S. (2005). Tropical atlantic influence on european heat waves. J. Clim. 18, 2805–2811. doi:10.1175/JCLI3506.1

CrossRef Full Text | Google Scholar

Dell’Aquila, A., Lucarini, V., Ruti, P. M., and Calmanti, S. (2005). Hayashi spectra of the northern hemisphere mid-latitude atmospheric variability in the NCEP–NCAR and ECMWF reanalyses. Clim. Dyn. 25, 639–652. doi:10.1007/s00382-005-0048-x

CrossRef Full Text | Google Scholar

Duchez, A., Frajka-Williams, E., Josey, S. A., Evans, D. G., Grist, J. P., Marsh, R., et al. (2016). Drivers of exceptionally cold north atlantic ocean temperatures and their link to the 2015 european heat wave. Environ. Res. Lett. 11, 074004. doi:10.1088/1748-9326/11/7/074004

CrossRef Full Text | Google Scholar

Fraedrich, K., Jansen, H., Kirk, E., Luksch, U., and Lunkeit, F. (2005). The planet simulator: towards a user friendly model. Meteorol. Z. 14, 299–304. doi:10.1127/0941-2948/2005/0043

CrossRef Full Text | Google Scholar

Fraedrich, K., Kirk, E., and Lunkeit, F. (1998). Puma: portable university model of the atmosphere. Dtsch. Klimarechenzentrum 38.

Google Scholar

Gálfi, V. M., and Lucarini, V. (2021). Fingerprinting heatwaves and cold spells and assessing their response to climate change using large deviation theory. Phys. Rev. Lett. 127, 058701. doi:10.1103/PhysRevLett.127.058701

PubMed Abstract | CrossRef Full Text | Google Scholar

Gálfi, V. M., Lucarini, V., Ragone, F., and Wouters, J. (2021). Applications of large deviation theory in geophysical fluid dynamics and climate science. La Riv. del Nuovo Cimento 44, 291–363. doi:10.1007/s40766-021-00020-z

CrossRef Full Text | Google Scholar

Gálfi, V. M., Lucarini, V., and Wouters, J. (2019). A large deviation theory-based analysis of heat waves and cold spells in a simplified model of the general circulation of the atmosphere. J. Stat. Mech. Theory Exp. 2019, 033404. doi:10.1088/1742-5468/ab02e8

CrossRef Full Text | Google Scholar

García-Herrera, R., Díaz, J., Trigo, R. M., Luterbacher, J., and Fischer, E. M. (2010). A review of the european summer heat wave of 2003. Crit. Rev. Environ. Sci. Technol. 40, 267–306. doi:10.1080/10643380802238137

CrossRef Full Text | Google Scholar

Gardiner, C. W. (1985). Handbook of stochastic methods, 3. springer Berlin.

Google Scholar

Haarsma, R. J., Selten, F. M., and Drijfhout, S. S. (2015). Decelerating atlantic meridional overturning circulation main cause of future west european summer atmospheric circulation changes. Environ. Res. Lett. 10, 094007. doi:10.1088/1748-9326/10/9/094007

CrossRef Full Text | Google Scholar

Hans von Storch, F. W. Z. (2002). Statistical analysis in climate research. Cambridge University Press.

Google Scholar

Hayashi, Y. (1971). A generalized method of resolving disturbances into progressive and retrogressive waves by space fourier and time cross-spectral analyses. J. Meteorological Soc. Jpn. Ser. II 49, 125–128. doi:10.2151/jmsj1965.49.2_125

CrossRef Full Text | Google Scholar

Hayashi, Y. (1979). A generalized method of resolving transient disturbances into standing and traveling waves by space-time spectral analysis. J. Atmos. Sci. 36, 1017–1029. doi:10.1175/1520-0469

CrossRef Full Text | Google Scholar

Herbert, C., and Bouchet, F. (2017). Predictability of escape for a stochastic saddle-node bifurcation: when rare events are typical. Phys. Rev. E 96, 030201. doi:10.1103/physreve.96.030201

PubMed Abstract | CrossRef Full Text | Google Scholar

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The era5 global reanalysis. Q. J. R. Meteorological Soc. 146, 1999–2049. doi:10.1002/qj.3803

CrossRef Full Text | Google Scholar

Horton, R. M., Mankin, J. S., Lesk, C., Coffel, E., and Raymond, C. (2016). A review of recent advances in research on extreme heat events. Curr. Clim. Change Rep. 2, 242–259. doi:10.1007/s40641-016-0042-x

CrossRef Full Text | Google Scholar

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., et al. (2013). The community earth system model: A framework for collaborative research. Bull. Am. Meteorological Soc. 94, 1339–1360. doi:10.1175/BAMS-D-12-00121.1

CrossRef Full Text | Google Scholar

IPCC (2014). Climate change 2013 – the physical science basis: Working group I contribution to the fifth assessment report of the intergovernmental panel on climate change. Cambridge University Press. doi:10.1017/CBO9781107415324

CrossRef Full Text | Google Scholar

Jacques-Dumas, V., Ragone, F., Borgnat, P., Abry, P., and Bouchet, F. (2022). Deep learning-based extreme heatwave forecast. Front. Clim. 4. doi:10.3389/fclim.2022.789641

CrossRef Full Text | Google Scholar

Kornhuber, K., Coumou, D., Vogel, E., Lesk, C., Donges, J. F., Lehmann, J., et al. (2020). Amplified rossby waves enhance risk of concurrent heatwaves in major breadbasket regions. Nat. Clim. Change 10, 48–53. doi:10.1038/s41558-019-0637-z

CrossRef Full Text | Google Scholar

Kornhuber, K., Osprey, S., Coumou, D., Petri, S., Petoukhov, V., Rahmstorf, S., et al. (2019). Extreme weather events in early summer 2018 connected by a recurrent hemispheric wave-7 pattern. Environ. Res. Lett. 14, 054002. doi:10.1088/1748-9326/ab13bf

CrossRef Full Text | Google Scholar

Lau, W. K. M., and Kim, K.-M. (2012). The 2010 Pakistan flood and Russian heat wave: teleconnection of hydrometeorological extremes. J. Hydrometeorol. 13, 392–403. doi:10.1175/JHM-D-11-016.1

CrossRef Full Text | Google Scholar

Lestang, T., Ragone, F., Bréhier, C.-E., Herbert, C., and Bouchet, F. (2018). Computing return times or return periods with rare event algorithms. J. Stat. Mech. Theory Exp. 2018, 043213. doi:10.1088/1742-5468/aab856

CrossRef Full Text | Google Scholar

Miloshevich, G., Cozian, B., Abry, P., Borgnat, P., and Bouchet, F. (2023). Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data. Phys. Rev. Fluids 8, 040501. doi:10.1103/PhysRevFluids.8.040501

CrossRef Full Text | Google Scholar

Muñoz Sabater, J., Dutra, E., Agustí-Panareda, A., Albergel, C., Arduini, G., Balsamo, G., et al. (2021). Era5-land: A state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data 13, 4349–4383. doi:10.5194/essd-13-4349-2021

CrossRef Full Text | Google Scholar

Osman, M. B., Coats, S., Das, S. B., McConnell, J. R., and Chellman, N. (2021). North atlantic jet stream projections in the context of the past 1,250 years. Proc. Natl. Acad. Sci. 118, e2104105118. doi:10.1073/pnas.2104105118

PubMed Abstract | CrossRef Full Text | Google Scholar

Otto, F., Massey, N., Van Oldenborgh, G. J., Jones, R., and Allen, M. (2012). Reconciling two approaches to attribution of the 2010 Russian heat wave. Geophys. Res. Lett. 39, 4702. doi:10.1029/2011GL050422

CrossRef Full Text | Google Scholar

Perkins, S. E. (2015). A review on the scientific understanding of heatwaves—Their measurement, driving mechanisms, and changes at the global scale. Atmos. Res. 164-165, 242–267. doi:10.1016/j.atmosres.2015.05.014

CrossRef Full Text | Google Scholar

Petoukhov, V., Petri, S., Kornhuber, K., Thonicke, K., Coumou, D., and Schellnhuber, H. J. (2018). Alberta wildfire 2016: apt contribution from anomalous planetary wave dynamics. Sci. Rep. 8, 12375. doi:10.1038/s41598-018-30812-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Petoukhov, V., Petri, S., Rahmstorf, S., Coumou, D., Kornhuber, K., and Schellnhuber, H. J. (2016). Role of quasiresonant planetary wave dynamics in recent boreal spring-to-autumn extreme events. Proc. Natl. Acad. Sci. 113, 6862–6867. doi:10.1073/pnas.1606300113

PubMed Abstract | CrossRef Full Text | Google Scholar

Petoukhov, V., Rahmstorf, S., Petri, S., and Schellnhuber, H. J. (2013). Quasiresonant amplification of planetary waves and recent northern hemisphere weather extremes. Proc. Natl. Acad. Sci. 110, 5336–5341. doi:10.1073/pnas.1222000110

PubMed Abstract | CrossRef Full Text | Google Scholar

Philip, S., Kew, S., van Oldenborgh, G. J., Otto, F., Vautard, R., van der Wiel, K., et al. (2020). A protocol for probabilistic extreme event attribution analyses. Adv. Stat. Climatol. Meteorology Oceanogr. 6, 177–203. Publisher: Copernicus GmbH. doi:10.5194/ascmo-6-177-2020

CrossRef Full Text | Google Scholar

Philip, S. Y., Kew, S. F., van Oldenborgh, G. J., Anslow, F. S., Seneviratne, S. I., Vautard, R., et al. (2021). Rapid attribution analysis of the extraordinary heatwave on the pacific coast of the us and Canada june 2021. Earth Syst. Dyn. Discuss. 2021, 1–34. doi:10.5194/esd-2021-90

CrossRef Full Text | Google Scholar

Pratt, R. W. (1976). The interpretation of space-time spectral quantities. J. Atmos. Sci. 33, 1060–1066. doi:10.1175/1520-0469(1976)033⟨1060:TIOSTS⟩2.0.CO;2

CrossRef Full Text | Google Scholar

Quemener, E. (2014). SIDUS-the solution for extreme deduplication of an operating system. Linux J. 2013.

Google Scholar

Ragone, F., and Bouchet, F. (2020). Computation of extreme values of time averaged observables in climate models with large deviation techniques. J. Stat. Phys. 179, 1637–1665. doi:10.1007/s10955-019-02429-7

CrossRef Full Text | Google Scholar

Ragone, F., and Bouchet, F. (2021). Rare event algorithm study of extreme warm summers and heatwaves over europe. Geophys. Res. Lett. 48, e2020GL091197. doi:10.1029/2020gl091197

CrossRef Full Text | Google Scholar

Ragone, F., Wouters, J., and Bouchet, F. (2018). Computation of extreme heat waves in climate models using a large deviation algorithm. Proc. Natl. Acad. Sci. 115, 24–29. doi:10.1073/pnas.1712645115

PubMed Abstract | CrossRef Full Text | Google Scholar

Rousi, E., Kornhuber, K., Beobide-Arsuaga, G., Luo, F., and Coumou, D. (2022). Accelerated western european heatwave trends linked to more-persistent double jets over eurasia. Nat. Commun. 13, 3851. doi:10.1038/s41467-022-31432-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Rousi, E., Selten, F., Rahmstorf, S., and Coumou, D. (2021). Changes in north atlantic atmospheric circulation in a warmer climate favor winter flooding and summer drought over europe. J. Clim. 34, 2277–2295. doi:10.1175/JCLI-D-20-0311.1

CrossRef Full Text | Google Scholar

Schiermeier, Q. (2021). Climate change made Europe's mega-heatwave five times more likely. Nature 571, 155. Bandiera_abtest: a Cg_type: News Publisher: Nature Publishing Group Subject_term: Climate change, Climate sciences. doi:10.1038/d41586-019-02071-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Schubert, S., Wang, H., and Suarez, M. (2011). Warm season subseasonal variability and climate extremes in the northern hemisphere: the role of stationary rossby waves. J. Clim. 24, 4773–4792. doi:10.1175/JCLI-D-10-05035.1

CrossRef Full Text | Google Scholar

Seneviratne, S., Zhang, X., Adnan, M., Badi, W., Dereczynski, C., Di Luca, A., et al. (2021). Weather and climate extreme events in a changing climate. Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press, 1513–1766. doi:10.1017/9781009157896.013

CrossRef Full Text | Google Scholar

Stéfanon, M., D’Andrea, F., and Drobinski, P. (2012). Heatwave classification over europe and the mediterranean region. Environ. Res. Lett. 7, 014023. doi:10.1088/1748-9326/7/1/014023

CrossRef Full Text | Google Scholar

Teng, H., Branstator, G., Wang, H., Meehl, G., and Washington, W. (2013). Probability of us heat waves affected by a subseasonal planetary wave pattern. Nat. Geosci. 6, 1056–1061. doi:10.1038/ngeo1988

CrossRef Full Text | Google Scholar

Van Oldenborgh, G. J., Wehner, M. F., Vautard, R., Otto, F. E. L., Seneviratne, S. I., Stott, P. A., et al. (2022). Attributing and projecting heatwaves is hard: we can do better. Earth’s Future 10, e2021EF002271. doi:10.1029/2021EF002271

CrossRef Full Text | Google Scholar

World Weather Attribution (2021). Western north american extreme heat virtually impossible without human-caused climate change.

Google Scholar

Zhang, R., Sun, C., Zhu, J., Zhang, R., and Li, W. (2020). Increased european heat waves in recent decades in response to shrinking arctic sea ice and eurasian snow cover. npj Clim. Atmos. Sci. 3, 7. doi:10.1038/s41612-020-0110-8

CrossRef Full Text | Google Scholar

Zscheischler, J., Martius, O., Westra, S., Bevacqua, E., Raymond, C., Horton, R. M., et al. (2020). A typology of compound weather and climate events. Nat. Rev. Earth Environ. 1, 333–347. doi:10.1038/s43017-020-0060-z

CrossRef Full Text | Google Scholar