2.1 Atmospheric circulation

We use the reanalysis data of the National Centers for Environmental Prediction (NCEP) (Kistler et al., 2001). We consider the geopotential height at 500 mb (Z500) over the North Atlantic region for the computation of circulation analogues. Sea-level pressure (SLP) is used for a posteriori diagnostics. We used the daily averages between 1 January 1948 and 31 December 2018. The horizontal resolution is 2.5^∘ in longitude and latitude. The rationale of using this reanalysis is that it covers 70 years and is regularly updated.

We consider Z500/SLP fields over two regions outlined in Fig. 1. The region in red is used to compute Z500 analogues. It is similar to the one advocated by Jézéquel et al. (2018). The region in blue is used for verification. The reason to use Z500 for computations, rather than SLP, is linked to the “heat low” of SLP during heat waves (Jézéquel et al., 2018). SLP and Z500 yield similar properties during the winter.

One of the caveats of this reanalysis dataset is the lack of homogeneity of assimilated data, in particular before the satellite era. This can lead to breaks in pressure related variables, although such breaks are mostly detected in the Southern Hemisphere and the Arctic regions (Sturaro, 2003), and it marginally impacts the eastern North Atlantic region.

Figure 1Location of five European stations. Longitudes are expressed in degrees east. Latitudes are expressed in degrees north. The red rectangle indicate the zone on which Z500 circulation analogues are computed. The blue rectangle indicates the region of pattern verification for Z500 and SLP.

Since Z500 values depend on temperature, we detrend the Z500 daily field by removing a seasonal average linear trend from each grid point. This preprocessing is performed to ensure that the results do not depend on atmospheric trends. All the analogue computations of this paper were performed on detrended and raw Z500 data, in order to verify the robustness of the results to Z500 trends.

2.2 Temperature observations

We took daily averages (TGs) of temperatures from the European Climate Assessment & Dataset (ECA & D) project (Klein-Tank et al., 2002). We extracted data from Berlin, De Bilt, Toulouse, Orly and Madrid (Fig. 1). We consider data from June–July–August (JJA). Those five stations cover a large longitudinal and latitudinal range in western Europe. These datasets were also chosen because

• they start before 1948 and end after 2018. This allows the computation of analogue temperatures with the Z500 from the NCEP reanalysis, which includes that period;

• they contain less than 10 % of missing data.

These two criteria allow keeping 528 of the 11 422 ECA&D stations that are available in 2018.

The five daily TGs are then averaged in order to provide a daily European temperature index. This choice is made to simplify the presentation of results. Results for individual stations are presented in the Supplement. This would overcome the potential caveat of the European temperature index, which does not capture the spatial variability of European heat waves (Stefanon et al., 2012). Since we focus on extremely hot temperature spells, this proof-of-concept study is limited to summer, and we consider the period between 1 June to 31 August.

3.1 Circulation analogues

Analogues of circulation are computed on Z500 data from reanalysis data of the NCEP reanalysis (Sect. 2.1) on the region outlined in red in Fig. 1. The reason to use this field rather than SLP is linked to the heat low of SLP during heat waves (Jézéquel et al., 2018). SLP and Z500 yield similar properties during the winter. For each day between 1 January 1948 and 31 December 2018, the best 20 analogues (with respect to a Euclidean distance) in a different year and within 30 calendar days are searched. This follows the procedure of Yiou et al. (2013).

The small region (20^∘ W–25^∘ E; 40–70^∘ N) is used to compute analogues of Z500, to simulate continental temperatures, following the domain recommendations of the analysis of Jézéquel et al. (2018). We compute analogues on “raw” daily Z500 data and detrended Z500. Detrending is performed by removing a smoothing spline on the spatially averaged Z500. The increasing trend of surface temperatures necessitates this trend removal in Z500.

The larger region (80^∘ W–25^∘ E; 30–70^∘ N) covers the North Atlantic region. The atmospheric circulation of the North Atlantic evolves within this region. This region is used to compute composites of large-scale SLP/Z500 fields during simulated events, for an evaluation of large-scale features during those small-scale events.

The analogues of Z500 are computed with the “blackswan” Web Processing Service (WPS) described by Hempelmann et al. (2018).

3.2 Stochastic weather generators

Ensembles of simulations of temperature can be performed with the rules illustrated by Yiou (2014), with analogue-based stochastic weather generators (SWGs). This type of SWG is equivalent to a resampling procedure (Ailliot et al., 2015).

A static SWG can be defined to simulate surrogate ensembles of Z500 sequences that are analogous to observed Z500 sequences. For each day t between 1 June and 31 August, we keep the K=20 best Z500 analogues. We randomly select one day (k) among those K analogues and t (i.e., among K+1 days), with a probability weight that is inversely proportional to the correlation of the K analogues with the Z500 pattern at time t. This constraint favors analogues with the best patterns among those with the closest distance. This also favors the choice of t. With this type of SWG, simulated trajectories are random perturbations (by analogues) of an observed trajectory.

A dynamic SWG is defined to simulate ensembles of Z500 sequences that could have been, from a given initial condition. For an initial day t (e.g., a 1 June), we have K best Z500 analogues. We randomly select one date $\stackrel{\mathrm{̃}}{t}$ among the dates of the K analogues and t (hence K+1 dates), with a probability weight that is

1. inversely proportional to the number of calendar days between the analogues dates $\stackrel{\mathrm{̃}}{t}$ and t. This constrains the time of analogues to move forward. The weights of the analogues can be chosen proportional to $\exp \left(-{\mathit{\alpha }}_{\mathrm{cal}}|\stackrel{\mathrm{̃}}{t}-t|\right)$, where α_cal≥0 weighs the importance given to seasonality or the calendar day;

2. proportional to the correlation of the analogue with the Z500 pattern at time t. This constraint favors analogues with the best patterns among those with the closest distance.

The simulated next day t^′ (e.g., a 2 June) is then the next day of the selected analogue: $t^{\prime }=\stackrel{\mathrm{̃}}{t}+\mathrm{1}$ of the initial day (1 June). Then t is replaced by t^′. This random selection of analogues is sequentially repeated until a lead time T=90 d to simulate a whole summer. This generates one random daily trajectory of Z500 or any climate variable between t and t+T. Note that the dynamic weather generator spans a wider range of possibilities than the static weather generator, which is constrained by observed trajectories.

These two types of SWG (static and dynamic) give useful information that is exploited in this paper. The random sampling procedures are repeated S times to generate an ensemble of trajectories. Yiou (2014) showed that temperature biases are rather small but the time auto-covariance is underestimated in both SWGs.

3.3 Empirical importance sampling

The idea behind importance sampling is to simulate trajectories of a physical system that optimize a criterion in a computationally efficient way. Ragone et al. (2017) used such an algorithm to simulate extreme heat waves with an intermediate complexity climate model. The procedure of importance sampling algorithms, say to simulate extreme heat waves with a climate model, is to start from an ensemble of S initial conditions and compute trajectories of the climate model from those initial conditions. An optimization observable is defined for the system. In this case, it can be the spatially averaged temperature over Europe. The trajectories for which the observable (daily average temperature) is lowest during the first steps of simulation are deleted and replaced by small perturbations of remaining ones. In this way, each time increment of the simulations keeps trajectories with the highest values of the observable. At the end of the season, one obtains S simulations for which the observable (here average temperature over Europe) has been maximized. Since those trajectories are solutions of the equations of a climate model, they are necessarily physically consistent. Ragone et al. (2017) argue that the probability of the simulated trajectories is controlled by a parameter that weighs the importance to the highest observable values: if n trajectories are deleted at each time step, the simulation of an ensemble of T-long trajectories has a probability of ${\left(\mathrm{1}-n/S\right)}^T$. Hence one obtains a set of S trajectories with very low probability after T time increments, at the cost of the computation of S trajectories. For comparison purpose, if one wants to obtain S trajectories that have a low probability (p) observable, then the number of necessary “unconstrained” simulations is of the order of S∕p, so that most of those simulations are left out. Systems like weather@home (Massey et al., 2015) that generate tens of thousands of climate simulations are just sufficient to obtain S=100 centennial heat waves, and the number of “wasted” simulations is very high. Therefore, importance sampling algorithms are very efficient ways to circumvent this difficulty.

Figure 2Illustration of the analogue importance sampling for mean daily temperature (TG). The red rectangles indicate the selected Z500 analogues for date t⁽ⁱ⁾. The red stars in the lower graphs are the temperature values during the selected analogue date. Blue arrows indicate the correspondence between Z500 analogues and daily TG. Here, TG is the daily mean temperature for the Berlin–De Bilt–Orly–Toulouse–Madrid (BDOTM) average. Panel (b) illustrates simulated TG (red, orange and brown lines); the black line is the observed TG; the blue line is the seasonal cycle of TG.

Download

Here, we propose an adaptation of such an algorithm to the stochastic weather generators of Sect. 3.2 to simulate extreme heat waves. The observable to be optimized is average daily temperature (named TG in the ECA & D nomenclature, Klein-Tank et al., 2002). To this end, we propose a new “rule” to the SWGs (static and dynamic) by giving more weight to the hottest temperatures. At each time step t, we note t^(k) the dates of the K best analogues. The temperature values of the K analogues of t and TG_t are sorted in decreasing order. The ranks are written R_k ($k\in \mathit{\{}\mathrm{0},\mathrm{\dots },K\mathit{\}}$). For example, the rank of the hottest temperature among analogues and temperature at day t is 1. We chose weights ${\left\{w^{\left(k\right)}\right\}}_{k\in \mathit{\{}\mathrm{0},\mathrm{\dots },K\mathit{\}}}$ so that

where α is a positive number and A is a normalizing constant so that the sum of weights over k is 1:

The useful property of this formulation of weights is that the values of w^(k) do not depend on time t because the rank values R_k are integers between 1 and K+1. The weight values do not depend on the unit of the variable either, so that this procedure does not need major adaptation to simulate other types of climate variables (e.g., precipitation or wind speed). If α=0, this is equivalent to a stochastic weather generator described in Sect. 3.2.

The date for the next day is chosen at random by sampling $\mathit{\{}t,t^{\left(\mathrm{1}\right)},\mathrm{\dots },t^{\left(K\right)}\mathit{\}}$ with the weights ${\left\{w^{\left(k\right)}\right\}}_{k\in \mathit{\{}\mathrm{0},\mathrm{\dots },K\mathit{\}}}$. Therefore, if $\stackrel{\mathrm{̃}}{T}$ is the simulated temperature among $\mathit{T}=\mathit{\{}{\mathrm{TG}}_t,{\mathrm{TG}}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathrm{TG}}^{\left(K\right)}\mathit{\}}$, then $\mathrm{Pr}\left(\stackrel{\mathrm{̃}}{T}={\mathit{T}}^{\left(k\right)}\right)=w^{\left(k\right)}$. The expected value $E\left(\stackrel{\mathrm{̃}}{T}\right)$ is then

where sort(T) are the sorted values of T in descending order. This allows us to select the circulation analogues that favor the highest temperature. The α parameter gives some flexibility to select analogues with lower temperatures. It plays the same role as the simulation parameter of Ragone et al. (2017), which controls the return times of trajectories. This choice of weights is also interesting because it allows for an approximation of the mean value of simulations as a function of the parameter α, as one can recognize a discrete Laplace transform of the distribution of T in Eq. (3).

This new rule replaces the analogue correlation weights of the dynamic and static SWGs defined in Sect. 3.2. The calendar day weight has to be maintained in order to keep a seasonality, especially for the dynamic scheme. Hence, this allows simulating

1. the hottest season that could have been with a similar atmospheric circulation, with a static analogue SWG;

2. the hottest season that could have been in the same climate, with a dynamic analogue SWG.

The algorithm is illustrated in Fig. 2.

The dynamic SWG option is closest to the importance sampling experiments performed by Ragone et al. (2017). The major algorithmic difference is that importance sampling “eliminates” trajectories with “poor” properties of the observable, while the analogue importance sampling “favors” the trajectories with the “best” properties.

The analogue importance sampling uses a self-coherent dataset as a basis (here: observations for temperature and a reanalysis for Z500 analogues), ensuring that the simulated trajectories bear physical consistence between temperature and the atmospheric circulation. The physical relevance also requires that simulations follow a seasonal cycle. If one is interested in simulating a perpetual summer, then the weight value α_cal on the calendar day is not very important. For physically realistic simulations, the value of α_cal can be chosen to be the smallest for which the temperature simulation median yields a seasonality. This prevents the system from staying in a perpetual mid-August.

The calendar “nudging” parameter α_cal needs to be chosen with care so that the additional rule does not create un-physical simulations (i.e., with time going backwards or simulating a perpetual summer). The value of α_cal was estimated by trial and error, by taking the smallest value for which most (e.g., more than 70 %) of the dynamic simulations end with dates after the second half of August. In our case (summer temperature simulations), a parameter value of 5 is deemed reasonable for summer temperature simulations. The value of α_cal could be different if another warm season has to be simulated because the dynamics of temperature variations depends on the season.

Figure 3Percentage of dynamic simulations (from 100 simulations) of extreme summers for which the last day falls after 15 August, as a function of the parameter α_cal, between 1 and 10. The vertical dashed line is for α_cal=5; the horizontal dashed line is for 70 %.

Download

From the numerical point of view, if α=0.1, 15 to 18 analogues have a probability larger than 0.1 of being selected. If seasonal trajectories (N=90 d) are simulated, this means that more than 15⁹⁰≈10¹⁰⁰ different trajectories of warm seasons are possible.

From Eq. (3), the probability distribution of the simulations is linked to the value of α. A formulation of the expected probability can be obtained heuristically. For example, let Q be the smallest number so that

where ϵ>0 is a small number (for example 1∕N, where N is the number of simulations that are needed to observe one simulated event) and A was given by Eq. (2). Then the probability of dynamic trajectories with parameter value α is close to (Q∕K)^M, where M is the average number of independent days during the simulated season (M≈18 for a season of 90 d) and K is the number of analogues (K=20). Such a heuristic formulation is close to what is obtained by Ragone et al. (2017). Yet this formulation suffers from numerical problems for large values of α (larger than 0.5) and small values of K (like 20). An alternate empirical estimation of the probability distribution of the average of the trajectories is to consider that the seasonal average of temperature closely follows a Gaussian distribution. Then an empirical estimate of the probability is obtained by comparing the quantiles of the observed distribution of average temperatures and the mean of simulated averages. With a chosen value of K=20, the heuristic and empirical approximations give similar estimates of the probability (or return period) for values of α≤0.5. We shall keep the empirical estimate of probabilities or return periods in the rest of the paper. Refined heuristic estimates are left for future (more statistical) work.

We emphasize that the methodology shown in this section (or in Ragone et al., 2017) is quite general and not specific to summer heat waves, as other types of extremes like cold spells (Cattiaux et al., 2010) or long episodes of precipitation (Schaller et al., 2016) can be envisioned.

4.1 Interannual trends

The distribution of temperature simulations is shown in Fig. 4, with parameter α values of 0, 0.2, 0.5 and 1. The static SWG distributions (blue box plots) yield variations that are closely correlated (r=0.75) with the observed values (black lines). The dynamic SWG distributions (red box plots) have lower-amplitude variations, although they are also correlated with observed values (r=0.58). Figure 4 shows that the SWG distribution means increase with the α parameter. When α=0, there is no importance sampling, so that both types of simulations fluctuate around the mean value of observations (Fig. 4a). In that case, the record value (in 2003) or second-highest value (in 2018) of temperature is not reached by the S=100 simulations. The record low temperature (in 1956) is also barely reached.

We observe no trend in the simulated trajectories with α>0 for all years. This is explained by a tendency to select recent analogues in the importance sampling because of the general temperature trend. The amplitude of the simulations is higher for dynamical simulations than for static simulations because the simulated trajectories are not constrained by the observed ones and could hence behave in a rather different way.

As soon as α>0, the distributions of simulated temperature exceed the observed values (Fig. 4b–d) because the SWG tends to select analogues for which temperature is higher than observed, by construction. When α>0, the dynamic trajectories generally yield higher values than the static trajectories because the dynamic SWG chooses the warmest analogues among all analogues. We observe that when α≥0.5, then the dynamic trajectories are all above the 2003 record. Roughly half of the static trajectories with α=1 stand above the 2003 record.

When α>0, the static SWG generates similar circulation patterns that lead to warmer summer temperatures (blue box plots). The 2003 record value in observations always yields the warmest JJA simulations. 2018 is the second-hottest summer in the observations, but it is no longer the second when α>0, with the static SWG, as a few years with heat waves with analogue atmospheric circulation could have exceeded the 2018 temperature value.

We observe an increase in average extreme summers with parameter α. Figure 5 summarizes the temperature probability distributions for all years, in the observations and simulations, for varying values of α. Here, return periods are estimated from a Gaussian approximation of the variability of average temperature (e.g., the black box plot in Fig. 5) for present-day conditions. This is justified by noting that we consider the probability distribution of temporal averages (92 d in JJA) of spatial averages (five stations). Those averages motivate a Gaussian approximation from the central limit theorem (von Storch and Zwiers, 2001).

Figure 5Box plots of simulated June–July–August (JJA) temperatures (TG in degrees Celsius) as a function of the parameter α for all years. The boxes of box plots indicate the median (q₅₀), lower (q₂₅) and upper (q₇₅) quartiles of the simulated values. The upper whiskers indicate $min[max(T),\mathrm{1.5}\times (q_{\mathrm{75}}-q_{\mathrm{25}}\left)\right]$. The lower whisker has a symmetric formulation. The points are the simulated values that are above or below the defined whiskers.The axis on the right indicates return times (in years), assuming a Gaussian distribution of JJA temperature averages with parameters estimated from the white box plot. The horizontal dotted line is the mean TG in 2003.

Download

Here we obtain a return time of the 2003 summer heat wave of around 350 years in present-day conditions (i.e., the last decades) (Fig. 5). The return time of this heat wave has been estimated to be as large as 10⁶ years (Schaer et al., 2004; Stott et al., 2004; Cattiaux and Ribes, 2018). Such results were based on extrapolating the probability density law of temperature observations since 1900. Therefore the underlying hypothesis to estimate the return period (or probability) is different from the estimate of Schaer et al. (2004). The SWG approach answers the question of how likely the occurrence of an event (in present-day conditions) is. Note that the event itself does not need to have occurred. Using long time series answers the question of what the frequency of the event is. In that case, the event must have occurred. Hence, our estimate is closer to a “hitting time” than a “return time”, according to the definition of Haydn et al. (2005).

For the rest of the paper, we chose a value of α=0.5, which corresponds to heat waves whose intensity is comparable to the 2003 record for the static SWG and higher for the dynamic SWG (Fig. 5). The simulated heat waves are more intense than the one of 2003, although the climate conditions are similar to present-day ones. Therefore, such heat waves are even rarer than the 2003 record value.

For comparison purposes, we simulated local heat waves in each of the five stations. Figures A1, B1 and C1 in the Appendix show the probability distribution of simulated summers (with α=0.5) for Berlin, Orly and Madrid.

Figure 6Observed and simulated daily temperature (TG) variations (for the Berlin–De Bilt–Orly–Toulouse–Madrid (BDOTM) average) for the four selected summers: (a) 2003, (b) 1956, (c) 1986 and (d) 2018. Temperatures are expressed in degrees Celsius. The simulation parameters are α_cal=5 and α=0.5. The blue lines indicate the average seasonal cycle. The black lines are the observed TG. The red lines are the simulations with the highest JJA mean TG. The dashed lines are the simulations with the 95th quantile of JJA mean TG. The brown lines are for the median simulations.

Download

4.2 Daily variations

The daily variations in dynamic simulations of average temperature for the 4 selected years (1956, 1986, 2003 and 2018) are shown in Fig. 6, with α=0.5. This figure illustrates how the simulated temperatures evolve above the seasonal cycle.

The values exceed the observed values within 3 d, especially for the coldest summer (in 1956, Fig. 6b). The maximum of the ensembles (red lines) generally reaches a plateau in 10 d after 1 June. The median and 95th quantile trajectories yield a higher temporal variability, which is lower than the observed one. We note that the algorithm allows selecting colder Z500 analogues, so that temperature does not continuously increase but yields time variability. This also means that all trajectories are different from each other.

The daily variations also explain the weak (visual) correlation between the means of simulations and observations in Fig. 4. The whole correlation is constrained by the few days after the initial condition (1 June), as the simulations progressively “forget” the initial condition, especially for the cold summer in 1956.

4.3 Atmospheric patterns

We focus on the atmospheric circulation patterns that prevail in summer (JJA) over the North Atlantic region for simulation ensembles, during the 4 selected years. We show the SLP (Fig. 7) and Z500 (Fig. 8) composites over JJA. The composites are the JJA averages for reanalyses and over all static and dynamic simulations.

The observed mean JJA SLP patterns are very different for each year (Fig. 7a–d). In particular, the two major heat waves of 2003 (hottest) and 2018 (second-hottest) shown in Fig. 4 yield rather contrasting SLP and Z500 pattern anomalies at the North Atlantic scale, although both years are characterized by positive SLP/Z500 anomalies over western Europe.

Static simulations show mean JJA SLP and Z500 patterns that are similar to the observed ones for 2003, 1956, 1986 and 2018 (Figs. 7e–h and 8e–h). The mean SLP/Z500 patterns show deepened or shifted structures to maximize mean temperature. Those figures illustrate how a slight modification of the atmospheric structure could increase European summer temperature. Therefore, small perturbations of the daily atmospheric Z500 structures (a few meters) can be associated with ≈4 K to the mean JJA temperature of the coldest summer in the time series (1956), with α=0.5. This temperature change is larger than the expected increase from a perfect gas law (≈0.4 K). Interestingly, the shifts in SLP patterns are much larger than for Z500.

Figure 7Maps of SLP anomaly composites (in Pa, with 50 Pa increments) for the warmest summer (2003), coolest (1956), median (1986) and 2018. The 0 hPa isoline is indicated with a thick blue line. Negative anomalies are in dashed isolines. Horizontal axes are for longitudes in degrees east. Vertical axes are for latitudes in degrees north. Upper row (a–d): mean SLP from NCEP reanalyses. Center row (e–h): static simulations. Bottom row (i–l): dynamic simulations.

Dynamic simulations simulate optimal atmospheric patterns that are very similar across the North Atlantic, with strong anticyclonic patterns over western Europe (Figs. 7i–l and 8i–l). The resulting SLP and Z500 patterns are very similar for the 4 selected years, with positive anomalies above western Europe. The Z500 conditions for 2018 contrast with the other years with positive anomalies in the central North Atlantic. This shows how the dynamic simulations differ from each other and depend on the initial conditions. The sensitivity to initial conditions was exploited by Yiou and Déandréis (2019) for ensemble forecasts with analogues.

The differences between the static and dynamic simulations illustrate the different concepts that those two weather generators convey. The dynamic simulations end up with rather similar ranges of temperatures and resembling SLP/Z500 structures over the whole North Atlantic (high pressure or Z500 over western Europe), which can be interpreted as an optimal pattern leading to major heat waves.

Figure 8Z500 anomaly composites (in meters, with 10 m increments) for the warmest summer (2003), coolest (1956), median (1986) and 2018, for European averages. The 0 m isoline is indicated with a thick blue line. Negative anomalies are in dashed isolines. Horizontal axes are for longitudes in degrees east. Upper row (a–d): mean Z500 from NCEP reanalyses. Center row (e–h): static simulations. Bottom row (i–l): dynamic simulations.

For comparison purposes, Figs. A2, B2 and C2 show the atmospheric Z500 patterns linked to local heat waves centered on Berlin, Orly and Madrid (from east to west). The optimal Z500 patterns with dynamic simulations for Orly are similar to the ones of the average temperature across the five stations.