# Maximal reachable temperatures for Western Europe in current climate

Robin Noyelle, Yi Zhang, Pascal Yiou, Davide Faranda

► To cite this version:

Robin Noyelle, Yi Zhang, Pascal Yiou, Davide Faranda. Maximal reachable temperatures for Western Europe in current climate. Environmental Research Letters, 2023, 10.1088/1748-9326/acf679. hal-04210152

HAL Id: hal-04210152

https://hal.science/hal-04210152v1

Submitted on 28 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
LETTER • OPEN ACCESS

# Maximal reachable temperatures for Western Europe in current climate

To cite this article: Robin Noyelle et al 2023 *Environ. Res. Lett.* 18 094061

View the [article online](#) for updates and enhancements.

## You may also like

- *Internal variability in European summer temperatures at 1.5 °C and 2 °C of global warming*
  Laura Suarez-Gutierrez, Chao Li, Wolfgang A Müller et al.

- *Avoided population exposure to extreme heat under two scenarios of global carbon neutrality by 2050 and 2060*
  Yadong Lei, Zhili Wang, Xiaoye Zhang et al.

- *Escalating global exposure to compound heat-humidity extremes with warming*
  Dawei Li, Jiacan Yuan and Robert E Kopp
ENVIRONMENTAL RESEARCH
LETTERS

<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

21 June 2023

REVISED

4 August 2023

ACCEPTED FOR PUBLICATION

4 September 2023

PUBLISHED

18 September 2023

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

LETTER

Maximal reachable temperatures for Western Europe in current climate

Robin Noyelle<fnref n="1" />,*, Yi Zhang<fnref n="2" />,<fnref n="3" />, Pascal Yiou<fnref n="1" /> and Davide Faranda<fnref n="1" />,<fnref n="4" />,<fnref n="5" />

E-mail: robin.noyelle@lsce.ipsl.fr

Keywords: heatwaves, climate change, maximum temperatures, flow analogues

Abstract

Human bodies, ecosystems and infrastructures display a non-linear sensibility to extreme temperatures occurring during heatwave events. Preparing for such events entails to know how high surface air temperatures can go. Here we examine the maximal reachable temperatures in Western Europe. Taking the July 2019 record-breaking heatwave as a case study and employing a flow analogues methodology, we find that temperatures exceeding 50 °C cannot be ruled out in most urban areas, even under current climate conditions. We analyze changes in the upper bound of surface air temperatures between the past (1940–1980) and present (1981–2021) periods. Our results show that the significant increase in daily maximum temperatures in the present period is only partially explained by the increase of the upper bound. Our results suggest that most of the warming of daily maximum surface temperatures result from strengthened diabatic surface fluxes rather than free troposphere warming.

**1. Introduction**

Extreme weather events have strong impacts on soci-
eties and ecosystems [1]. Among them, heatwaves
have been the focus of extensive attention due to
their increasing frequency and intensity with global
warming [2]. Most biological and technical systems
have a limiting capacity to cope with extreme tem-
peratures and can respond non linearly once cer-
tain thresholds are reached. These include human
bodies [3–6], plants [7, 8], ecosystems [9] and
infrastructures [10]. Urban areas are particularly at
risk due to the so-called heat island effect [11] which
enhances maximum temperatures during the day and
minimum temperatures during the night. The capa-
city to predict and anticipate future maximum intens-
ities of heatwaves is therefore of primary importance
for adaptation to climate change.

The physical mechanisms leading to heatwaves in
the midlatitudes combine specific atmospheric syn-
optic conditions and anomalously low soil moisture,
which can interact to further increase the intensity

of the event [11–14]. Heatwaves are associated with
a slow moving, sometimes called quasi-stationary,
high-amplitude Rossby wave [15] which is often
embedded in a hemispheric pattern of wave number
5 or 7 [16–18]. Above the heatwave region, a blocking
anticyclone builds up at mid- and upper-level tropo-
sphere in conjunction with a change in the jet stream's
climatological path, towards a large poleward meridi-
onal meander. The anticyclone sustains the poleward
advection of warm air along its western flank, adia-
batic warming by subsidence and clear skies at its cen-
ter. When linked to abnormally dry soils, these con-
ditions promote warming through short-wave insol-
ation, resulting in the allocation of incoming solar
energy towards sensible rather than latent heat [19–
26]. If the anticyclone stays stationary, extreme tem-
perature can be reached. Although it is well estab-
lished that global warming renders the temperat-
ure distribution hotter on average, how the coupled
dynamics between the atmosphere, the oceans and
the soils will evolve in response to this forcing during
the hottest days is still debated [11, 14, 27–31].

<footnote n="1">¹ Laboratoire des Sciences du Climat et de l'Environnement, UMR 8212 CEA-CNRS-UVSQ, Université Paris-Saclay & IPSL, Gif-sur-Yvette 91191, France</footnote>

<footnote n="2">² Department of Earth and Planetary Science, University of California, Berkeley, Berkeley, CA 94720, United States of America</footnote>

<footnote n="3">³ Miller Institute for Basic Research in Science, University of California, Berkeley, Berkeley, CA 94720, United States of America</footnote>

<footnote n="4">⁴ London Mathematical Laboratory, 8 Margravine Gardens London, W6 8RH London, United Kingdom</footnote>

<footnote n="5">5 Laboratoire de Météorologie Dynamique, Ecole Normale Supérieure, PSL Research University & IPSL, Paris, France</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
One way to quantify the risk associated with heatwaves is to measure the tail probabilities of the temperature distribution, e.g. the risk of the maximum annual temperature reaching a certain level. A classical approach to deal with this assessment is extreme value theory (EVT). EVT has been used to determine statistical models for maxima (or minima) of climate variables of interest [32], such as temperatures. It is based on a convergence principle of the probability distribution of block maxima or peaks-over-threshold [33–35] towards a generalized extreme value (GEV) distribution. It allows to compute return values corresponding to very large return periods (i.e. longer than the period of observations), even in non stationary contexts [36]. This approach is for example extensively used by attribution methods to compare the probabilities of reaching extreme temperatures between a counterfactual world without climate change and a factual world with climate change [37–39]. This allows to estimate how climate change made a particular event more (or less) likely.

Fitting a GEV distribution to extreme temperatures usually gives estimates of the shape parameter that are robustly negative (e.g. [40, 41]). This implies that, in a stationary context, the distribution of annual maximum temperatures is bounded upwards. This upper bound is often the information asked to the scientific community by practitioners who want to prepare for the worst case scenario. However, recent intense heatwaves such as in 2010 in Russia [18], 2019 in Western Europe [42] and 2021 in the Pacific Northwest [43, 44] have challenged the reliability of this statistical upper bound by breaking it sometimes by a large margin. Another approach to estimate maximal reachable temperatures would be to run a long simulation with a climate model to sample more extreme events than those observed in the past. Apart from the inherent limitations of models to represent correctly the entire temperature distribution—especially at its tail [40, 45, 46] — this method is limited by the fact that estimating precisely small probabilities requires an extremely long simulation [47], which is too costly for most complex models. Various approaches have been suggested to address this problem, like rare events algorithms [48–50] or ensemble boosting [51].

Recently, [52] proposed a physics-based theory to address the question of the upper bound for extreme temperatures over midlatitude continental areas. They provide evidence for the hypothesis that moist convective instability sets a fundamental upper limit for surface temperatures. This allows for an estimate of maximal temperatures based on large scale circulation of the free-troposphere and surface specific humidity. In this paper, we make use of this theory to investigate the maximal reachable temperature in Western Europe in the current climate. We use

the flow analogues attribution method [53–55]. We condition the estimate of maximal temperatures on the free-troposphere large scale circulation observed during the most intense event of the recent years in Western Europe, namely the July 2019 heatwave. 2019 saw two exceptional heatwaves in June and July [56] with maximum temperatures very unlikely in a world without anthropogenic climate change [45, 57].

We focus on the July 2019 heatwave—when most records of Western Europe cities were broken—and we estimate the maximal temperatures this event could have reached in the present climate. We also investigate how these upper bounds have changed with respect to the past and detail the physical mechanisms related to the changes in daily maximum of surface temperatures.

This paper is organized as follows: in section 2, we present the data used, the computation of the heatwave maximal temperature and the flow analogs methodology. Section 3 presents and discusses the results. Finally, the conclusions drawn from our analysis are presented in section 4.

## 2. Data and methods

We use the ERA5 reanalysis data set of the European Center for Medium Range Weather Forecasting (ECMWF) over the period 1940–2021 [58]. The data have a horizontal resolution of 0.25° × 0.25°. For illustration purposes, we additionally use daily temperature data from the Paris Montsouris weather station taken from the Infoclimat website (https://www.infoclimat.fr/climatologie/globale/paris-montsouris/07156.html).

To study how more intense an observed heatwave event could have been, we use the criteria developed by [52] who provided evidence for the hypothesis that convective instability fundamentally limits maximal surface air temperatures over midlatitude land. This hypothesis states that, during heatwaves, surface air temperature $T$ reaches its maximal value when the air column is neutrally stratified. Accordingly, the heatwave stops when precipitation is triggered by the convective instability of the boundary layer. This happens when the surface air parcel moist static energy

$$MSE = c_p T + L_v q + gz_s \quad (1)$$

is equal to the free-tropospheric saturation moist static energy $MSE_a^*$. Here $c_p$ is the specific heat of air at constant pressure, $L_v$ the latent heat of vaporization, $q$ the specific humidity at the surface, $g$ the gravitational acceleration and $z_s$ the geopotential height of the surface. We choose to define the free-tropospheric situation as the state of the atmosphere at 500 hPa and $MSE_a^*$ is therefore computed from equation (1) by
replacing *T* with the air temperature at 500 hPa *T*<sub>500</sub>, *q* with the saturation specific humidity at *T*<sub>500</sub>:

$$q_{\text{sat}}(T_{500}) \approx \frac{\epsilon e_{\text{sat}}(T_{500})}{500 \text{ hPa}} \quad (2)$$

where *ϵ* is the molar ratio between water vapor and dry air and *e*<sub>sat</sub>(*T*<sub>500</sub>) is the saturation vapor pressure given by the Clausius–Clapeyron relation—and *z*<sub>s</sub> with the geopotential height at 500 hPa *z*<sub>500</sub>. The maximal surface temperature is therefore:

$$T_{s,\max} = T_{500} + \frac{L_v}{c_p}(q_{\text{sat}}(T_{500}) - q) + \frac{g}{c_p}(z_{500} - z_s). \quad (3)$$

This temperature is maximized if *q* is null. However, such a low specific humidity is unrealistic for midlatitude land because of the presence of surface waters and vegetation. In the following, except specified explicitly, we therefore use the hourly June–July–August–September (JJAS) minimum *q*<sub>min</sub> ever observed over the period 1940–2021 at each grid point. Using *q* at 2 m to compute *T*<sub>s,max</sub> is a conservative choice because when an air parcel rises, it entrains surrounding dryer air, which decreases its specific humidity. For *T*<sub>500</sub> and *z*<sub>500</sub> we use the daily mean values as the diurnal cycle does not influence strongly these quantities.

*T*<sub>s,max</sub> defines a heatwave potential in the sense that it gives the maximal temperature potentially reachable given the synoptic situation. Whether this temperature will be reached is however not certain because, for example, the circulation pattern may change too fast for the temperature to have time to increase under the forcing of incoming energy fluxes. Letting *TX1d* be the observed daily temperature maximum at the surface, one can therefore define an index *I*<sub>TX1d</sub> for how much an observed heatwave is close to its maximal potential:

$$I_{\text{TX1d}} := \frac{\text{TX1d} - T_{\text{clim}}}{T_{s,\max} - T_{\text{clim}}}, \quad (4)$$

where *T*<sub>clim</sub> is a climatological temperature and *TX1d* is the observed daily maximum temperature. Here we use the July averaged daily mean temperature as a climatological value. When this index is close to 1, the event is close to its maximal potential intensity: it could not have been more intense.

To determine how global warming may have impacted the physical mechanisms leading to the extremely intense event of July 2019 in Western Europe, we employ an attribution methodology. The most common attribution method of a specific event [37] consists in comparing the probabilities of an observable reaching a certain level in a counterfactual world and in a factual world. The difference between the two periods usually lies in their global mean

temperature, i.e. worlds with or without anthropogenic global warming. The ratio of the probabilities gives how much more likely an observed event has been rendered because of global warming. This method however does not condition the result on the large scale circulation and therefore bears the risk of comparing atmospheric dynamics that are different even though they give rise to the same values for the observable of interest.

Here we make use of the flow analogues attribution methodology proposed by [55] and adapted from [53] (see also [59, 60]). The idea of the method is to compare the expectation of a variable of interest *X* in the counterfactual world (*F* = 0) and in the factual world (*F* = 1), conditional on the large scale circulation *C*(*ζ*) of the observed event *ζ*:

$$\Delta_{C(\zeta)}X = \mathbb{E}[X | F = 1, C(\zeta)] - \mathbb{E}[X | F = 0, C(\zeta)]. \quad (5)$$

In our analysis, *ζ* corresponds to the 25 July 2019 situation, the day where the heat peak was reached. The conditioning on the large scale circulation *C*(*ζ*) should be understood as being in a reasonable vicinity of the circulation of the event rather than as a strict equality, which cannot be enforced. This framework allows to answer the question: how a similar large scale circulation pattern in the two worlds leads to different outcomes in an observable of interest? If the difference *Δ*<sub>*C*</sub>*X* is statistically significant, then one can say that in the factual world the event has been rendered more (or less) intense by *Δ*<sub>*C*</sub>*X*. This method allows to disentangle the role of thermodynamics and dynamics (see [61] for the 2019 event). Whether the event has become more likely can in principle be obtained [53] if one can estimate the probabilities P[C(*ζ*) | *F*] and P[C(*ζ*) | *X* > *x*], where *x* is a given threshold, in the two worlds.

This method can be used with climate models output with or without anthropogenic forcing to define the factual and counterfactual worlds. However, models have known deficiencies, including biases [62] and incorrect dynamics of extremes under forcing over Western Europe [40, 45, 46]. As in [55] we therefore choose to rely on reanalysis outputs only, which, even though they also present limitations, are the closest we have to the actual past state of the atmosphere. Here, we define the counterfactual world as the past period (1940–1980) and the factual world as the present period (1981–2021).

To condition on the synoptic circulation, we select in both periods the best 40 analogs as the daily mean 500 hPa geopotential height (*z*<sub>500</sub>) maps minimizing the pointwise Euclidean distance with the daily mean *z*<sub>500</sub> map of the target day *ζ*. In the factual period we also impose that no analog can be found within a 20 calendar days window around the
25 July 2019 to avoid selecting the event itself. The reason for choosing to condition on $z_{500}$ is that in the quasi-geostrophic approximation relevant for synoptic scales, $z_{500}$ acts as an approximate streamfunction for the general circulation and has been shown to be a good conditioning field for temperatures [54]. Over the period 1940–2021, the troposphere has warmed and geopotential heights may therefore have inflated. To avoid this factor to bias the analogs selection, we remove a spatially uniform trend to the geopotential heights before computing the Euclidean distance. This ensures that the spatial gradients are preserved and therefore that analogs have a similar synoptic circulation. To compute the trend we average spatially the $z_{500}$ maps and fit a second degree polynomial to this time series. The second order was chosen because of the so-called warming pause in the global temperature of the Earth in the period 1950–1970 [63]. The analogs are computed over the domain depicted with a dashed box in figure 1 panel (c). This domain is chosen because it encompasses the observed heatwave region at the surface and the relevant associated synoptic pattern, namely a large poleward meander of the jet stream associated with a strong anticyclonic anomaly over Western Europe. The analogs are chosen only in the months June, July, August and September to ensure closeness with the event. We verified that the results are not qualitatively sensitive to the choice of the number of analogs provided that we extract 20 or 30 analogs maps.

Forty-one years long factual and counterfactual periods ensure a large enough statistical sample of $z_{500}$ maps to be able to select sufficiently representative analogs over such a limited geographical area. The counterfactual period can be assumed to satisfy the assumption of climate stationarity—for example with regards to the mean (not shown) or extreme (see e.g. figure 1 panel (a)) temperatures in summer—with respect to anthropogenic climate change. This period is also long enough for the interannual and interdecadal natural variability of the atmospheric circulation to be averaged over. The same is not necessarily true for the factual period for which there is a strong warming of Western Europe between the 1980s and the 2010s decades. A direct identification of the factual period and a stationary warmer world should therefore be made with care and it is more relevant to frame our results as showing a difference between a past world and a present world [64].

Once obtained the two sets of analogs for the factual and counterfactual periods, we average them out to search for significant differences in observables $X$ of interest. To determine significant changes, we apply a two-tail Welch t-test [65], with different variances, at each grid point. We mark as significant only grid point changes for which the p-value of the test is below 0.05. In the following, as the reference event takes place in July we chose to display anomalies with respect to the July monthly mean and

standard deviation over the period 1990–2021. Apart when specified explicitly, this convention concerns all anomalies.

## 3. Results and discussion

### 3.1. Description of the event and intensity potential

In the end of July 2019, Western Europe was hit by a heatwave which broke many records in France, Belgium, the Netherlands and Western Germany ([45], figure 1 panel (b)). The heatwave peaked on the 25 July and, for instance, the Paris Montsouris station recorded its largest temperature since 1900 with a new record of 42.6 °C, breaking the previous 1940 record by a margin of 2.2 °C (figure 1 panel (a)). The synoptic situation was characterized by a strongly meridional meander of the high level jet (figure 1 panel (c)) sustaining advection of warm air and anomalies of up to 7 K at 500 hPa above the North Sea (figure 1 panel (b)). Extreme surface temperatures were reached after an anomalously dry and hot month (30% of climatological precipitations and +2.4 °C temperature anomaly for the Paris region) following the previous June heatwave [56], which made the event develop under desiccated soils (figure 1 panel (d)). The event however only lasted a few days and daily maximum temperatures decreased by more than 10 K on the 26 July.

Figure 2 panel (a) shows the difference between the maximal surface air temperature $T_{s,max}$ and the observed daily maximum temperature TX1d for the 25 July 2019. To ensure realism in the value of $T_{s,max}$, we use the hourly JJAS minimum specific humidity $q_{min}$ ever observed over the period 1940–2021 at each grid point. Over the core heatwave region (Northern France and Benelux), the daily maximum temperature provided by the reanalysis data set is close to the maximal potential intensity by a margin of 4–6 K, and consequently the heatwave potential index is close to 90%.

Table 1 shows the estimation of $T_{s,max}$ for several major urban areas in Western Europe which all set their 1940–2021 maximum temperature record during the July 2019 event. For all of them, we used a four grid point average encompassing most of the main city and its surroundings. We compare the estimated maximal surface temperature to the maximum value of a non-stationary GEV distribution fitted on the summer temperature maxima over the 1940–2021 period using the ERA5 reanalysis data set (see supplementary materials for more information). We emphasize that these values are conservative with respect to what can be actually measured by a weather station at the surface because they are averaged over a roughly 50 km by 50 km square. For example, the value for TX1d for the Paris region with ERA5 gives 41.4 °C whereas the Paris Montsouris station measured 42.6 °C.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Annual maximum temperature for the Paris Montsouris station (black) and the grid point where the station is located for ERA5 (red). The black dashed line corresponds to the temperature reached on the 25 July 2019 at the Paris Montsouris station. (b) Daily maximum 2 m air temperature (colors) and daily mean 500 hPa temperature anomalies (contours). (c) Daily mean meridional wind at 250 hPa (colors) and geopotential height at 500 hPa (contours). (d) Daily mean normalized soil moisture anomaly (colors) and hourly mean 2 m specific humidity at the time when maximum 2 m air temperature is reached (hashed). The hashed regions represent anomalies below -0.5 standard deviation. All maps are drawn for the 25 July 2019.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Difference between the maximal air surface temperature $T_{s,max}$ and the observed daily maximum temperature $TX_{1d}$ for the 25 July 2019. (b) Heatwave potential index for the 25 July 2019. We show land grid points where the daily maximum temperature on the 25 July 2019 exceeded the 99th percentile of daily maximum temperature in July over the period 1990–2021. $T_{s,max}$ is computed using the hourly summer minimum specific humidity $q_{min}$ ever observed over the period 1941–2021 at each grid point.</figcaption></figure>

Strikingly, even in the current climate, it is not possible to rule out the possibility of reaching 50 °C in all of these cities once we include the uncertainty on the minimum surface specific humidity. For example, for the Paris area the ERA5 dataset provides a daily maximum of 41.4 °C, while $T_{s,max}$ is 46.6 °C [45.2, 53.8]. Reaching such temperatures would imply either a very intense drought, with very low values of surface specific humidity, and/or massive advection of heat from neighbouring regions. The non-stationary GEV estimate for the upper bound is significantly lower than $T_{s,max}$ by 5 to

10 K. The median estimate is even lower than the 2019 value for Paris and Lille. For most cities, the estimated $T_{s,max}$ is inside the uncertainty range of the GEV estimate, but much closer to the upper bound than to the lower one. This suggests that the statistical analysis of past time series of maximum temperatures underestimates the risk of reaching extreme temperatures.

## 3.2. Flow analogs analysis

Figure 3 panel (a) shows the difference between the factual and counterfactual periods of the analogs averaged daily maximum 2 m air temperature $TX_{1d}$
<table><caption>**Table 1.** Maximal temperatures in °C for several urban areas in Western Europe. All of them reached their 1940–2021 maximum temperature record during the July 2019 event. See supplementary materials for the detail of the computation of the uncertainty range.</caption><thead><tr><th>City</th><th>TX1d (ERA5)</th><th>T<sub>s,max</sub> (realistic)</th><th>T<sub>s,max</sub> (absolute)</th><th>GEV estimate</th></tr></thead><tbody><tr><td>Paris</td><td>41.4</td><td>46.6 [45.2, 54.6]</td><td>54.2 [53.8, 54.6]</td><td>40.8 [40.0, 47.2]</td></tr><tr><td>Lille</td><td>39.7</td><td>47.1 [45.0, 55.0]</td><td>54.9 [54.6, 55.0]</td><td>39.2 [38.3, 48.5]</td></tr><tr><td>Bruxelles</td><td>38.9</td><td>47.1 [45.4, 54.9]</td><td>54.8 [54.5, 54.9]</td><td>42.8 [37.6, 49.7]</td></tr><tr><td>Amsterdam</td><td>35.6</td><td>47.8 [44.7, 56.2]</td><td>55.9 [55.7, 56.2]</td><td>37.3 [33.7, 41.5]</td></tr><tr><td>Köln</td><td>40.2</td><td>48.4 [45.9, 55.5]</td><td>55.3 [55.0, 55.5]</td><td>42.1 [39.6, 48.4]</td></tr><tr><td>Franckfurt</td><td>39.3</td><td>50.1 [46.4, 56.8]</td><td>56.6 [56.3, 56.8]</td><td>42.2 [39.3, 49.2]</td></tr><tr><td>London</td><td>35.6</td><td>48.6 [46.9, 56.0]</td><td>55.8 [55.6, 56.0]</td><td>42.9 [36.3, 51.0]</td></tr><tr><td>Luxembourg</td><td>37.8</td><td>45.4 [45.2, 54.1]</td><td>53.8 [53.5, 54.1]</td><td>41.1 [37.6, 47.5]</td></tr></tbody></table>

<figure><img src="image_6.png" /></figure>

over Western Europe. In the heatwave region (France, England, Benelux and Western Germany) TX1d is between 2 and 3 K more intense in the recent period. It should also be noted that this analogs difference is stronger by approximately 1 K than the difference in the daily mean temperature (see figure C1 in supplementary materials).

Figure 3 panel (b) displays the difference of the absolute maximal surface temperature $T_{s,max}$. This quantity has not changed significantly over most of Western Europe—except a small region above the Netherlands which coincides with an increase of $T_{500}$ (not shown). Figure 3 panel (c) shows the percentage of change in TX1d explained by the increase in $T_{s,max}$. The latter corresponds to no more than 30% of the actual change between the two period over France (Paris region) but as much as 70% over Benelux and Western Germany.

The change in $T_{500}$ is the dominant factor for the stronger increase in absolute $T_{s,max}$ over the Netherlands, with an increase of 0.7 K, whereas only 0.4 K over Paris (see vertical profiles in figure C2). There is also a more significant drying of surface air over the Netherlands (figure C2 panels (b) and (d)). These two elements lead to a much stronger

increase in the analogs averaged lifting condensation level (LCL) between the two regions: 21 hPa over the Netherlands vs 14 hPa over Paris, although the difference between the two periods is non significant at the 5% level. Dryer air parcels in a hotter troposphere follow the dry adiabatic temperature profile longer, which leads to increased surface temperature as can be seen when comparing white and magenta lines in figure C2 panels (a) and (c). These elements suggest that under the synoptic atmospheric conditions associated with the 25 July 2019 event, most of the observed change in 2 m air temperature maximum in France between the two periods is caused by surface processes—either local diabatic heating or advection from neighboring regions.

Figure 4 shows the analogs analysis for TX1d and $z_{500}$ (first row), surface specific humidity $q$ and $T_{500}$ (second row). In both periods the averaged $z_{500}$ pattern matches the event $z_{500}$ pattern (figure 1 panel (c)) and the difference between the periods is minimal as can be seen by the fact that there are almost no significant differences in the $z_{500}$ field (figure 4 panel (c)). We are therefore confident on the quality of the analogs and the relevance of the comparison
<figure><img src="image_7.png" /><figcaption>**Figure 4.** First row: analogs averaged TX1d anomalies (colors) and daily mean $z_{500}$ (contours) for (a) the counterfactual and (b) the factual periods. (c) Difference between panels (a) and (b). For TX1d the anomalies are computed with respect to the 1990–2021 July monthly mean 2 m air temperature climatology. Second row: analogs averaged hourly mean 2 m specific humidity at the time when TX1d is reached (colors) and daily mean $T_{500}$ (contours) for (d) the counterfactual and (e) the factual periods. (f) Difference between panels (d) and (e). In the last column, for all fields only the grid point with a significant difference at the 5% level are displayed.</figcaption></figure>

between mean observables conditional on the synoptic circulation.

Figure 4 panel (f) shows no significant difference in the $T_{500}$ fields between the two periods. This should however be taken with care because in equation (3) this variable is used to compute $e_{sat}(T_{500})$ which evolves strongly non-linearly with $T_{500}$. A small difference in $T_{500}$, even non significant, can therefore lead to significant $T_{s,max}$ differences at the surface as shown in figure 3 panel (b). Finally, over the heatwave region the difference in specific humidity is not significant (figure 4 panel (f)), except over a small area in Western Germany.

In figure C3 in supplementary materials we present the same plot for upper layer soil moisture, which is a variable closely related to 2 m specific humidity (low soil moisture reduces the evaporation potential). There is a general drying pattern over most of the heatwave region in the factual period, but it is significant only over the same region as specific humidity in figure 4 panel (f).

The flow analogs method for these two variables is however less robust because soil moisture may show a much slower dynamics than $z_{500}$. The distribution of soil moisture (and therefore specific humidity) is likely to depend on the full trajectory followed by the atmospheric dynamics rather than its last state. The association between the synoptic circulation $C(\zeta)$ of the 25 July 2019 and dry soils/low surface specific humidity may therefore be a coincidence, especially for such a short extreme. As a consequence, it is possible that the framework presented in equation (5) may not be relevant for these variables.

Figure C4 in supplementary materials indeed shows that the JJAS daily minimum of hourly mean specific humidity at 2 m displays a significant

increase between the two periods for Western Europe (figure C4 panel (c)), probably as a result of moisture advection from nearby oceans [66]. This is apparently in contradiction with the analogs analysis presented here only if one dismisses the conditioning on the $z_{500}$ pattern. Nonetheless, whether the soil desiccation pattern observed between the two periods is a relevant feature of the change in the physical mechanisms leading to extreme temperatures in the two periods can be questioned.

Previous literature showed that there is a decrease of surface specific humidity during the hottest days [56, 67–69], but this may arise as both a cause and a consequence of intensifying surface temperature extremes [70]. We note also that although soil moisture is well represented in ERA5, the temporal drying trend may be underestimated [71].

Figure 5 shows the $T_{s,max}$ indicator computed in both periods using either the actual specific humidity at the surface $q$ or the 1940–2021 minimum specific humidity $q_{min}$. As already shown above, figure 5 panel (f) shows a significant change in absolute $T_{s,max}$ only over a limited region, which cannot explain the intensity change in TX1d between the two periods. Moreover, if one takes into account the moistening of minimum surface specific humidity over Western Europe in the factual period, this small difference may disappear completely as it would decrease $T_{s,max}$ in the factual period. When one uses the actual $q$ of the analogs, the difference is much stronger between the two periods: there is an increase of more than 3 K of $T_{s,max}$ over a large region encompassing most of Germany. This region however does not coincide perfectly with the region with a significant change in TX1d in figure 4 panel (c).
<figure><img src="image_8.png" /></figure>

**Figure 5.** First row: analogs averaged $T_{s,max}$ computed using the specific humidity of the analogs when TX1d is reached for (a) the counterfactual and (b) the factual periods. (c) Difference between panels (a) and (b). Second row: analogs averaged $T_{s,max}$ computed using the minimum specific humidity over the 1940–2021 period for (d) the counterfactual and (e) the factual periods. (f) Difference between panels (d) and (e). In the last column, the hashed grid points correspond to a significant difference at the 5% level.

<figure><img src="image_9.png" /></figure>

**Figure 6.** Analogues averaged of the sum of downwards energy fluxes at the surface for (a) the counterfactual and (b) the factual periods. (c) Difference between the two periods. In panel (c), the hashed grid points correspond to a significant difference at the 5% level.

To investigate this question, we analyze incoming energy fluxes. Even though the heatwave potential intensity has not changed much between the two periods, the actual capacity of a synoptic circulation to reach this maximal intensity may be greater in the factual period. Figure 6 displays the analogues average of the sum of daily mean incoming energy fluxes at the surface (downward positive). It includes latent and sensible heat fluxes, in addition to net shortwave radiation and downward longwave radiation. Figure 6 panel (c) shows that there is a significant and strong difference in incoming energy fluxes between the two periods which coincide with the heatwave region where the difference in TX1d is significant. These fluxes are strongly coupled to the boundary layer and soils states—especially with regards to the partitioning between sensible and latent heat. Therefore one

cannot claim that energy fluxes are a fully external cause of the increase of TX1d between the periods. However, they are an indicator of the strong increase in local diabatic heating.

The observed increase in TX1d could also be the result of increased sensible heat advection from neighboring regions with increasingly desiccated land surface [24, 72]. It is possible that this phenomenon of non-local self-propagation of dry and hot conditions played an important role for the 25 July 2019 event as shown by [56], especially with regards to the aftermath of the June heatwave which desiccated soils in the south of France and the Iberian peninsula. However, it is unlikely to be the case for all of the 40 analogues in both periods: as shown by [73], in Europe, hot events are generally associated with weaker horizontal transport, but strong adiabatic
## Authors' contributions

R N did the data collection and analysis. All authors contributed to the design of the analysis and the writing of the paper.

## Appendix A. Non-stationary GEV estimation

If we consider the maximum $M_n = \max\{X_1, \dots, X_n\}$ of n i.i.d. random variables $X_1, \dots, X_n$, the extremal theorems [33, 34] state that if there exist sequences $a_n > 0$ and $b_n$ so that the distribution $\mathbb{P}(\frac{M_n-b_n}{a_n} \le z)$ converges towards $G(z)$ when n goes to infinity, then $G(z)$ is given by:

$$G(z) = \begin{cases} \exp\left\{-\left[1 + \xi\left(\frac{z-\mu}{\sigma}\right)\right]^{-1/\xi}\right\} & \text{if } \xi \neq 0 \\ \exp\left\{-\exp\left[-\left(\frac{z-\mu}{\sigma}\right)\right]\right\} & \text{if } \xi = 0 \end{cases} \quad (A.1)$$

where $\mu$ is called the location parameter, $\sigma$ the scale parameter and $\xi$ the shape parameter. This distribution is usually called the GEV distribution. Importantly, it is bounded upwards for $\xi < 0$ and the upper bound is given by:

$$z_{\max} = \mu - \frac{\sigma}{\xi}. \quad (A.2)$$

Here we estimate a non-stationary GEV by letting the parameter $\mu$ depend on another covariate $T$ encompassing the warming of the Earth: $\mu = \mu_1 + \mu_2 T$. We chose $T$ to be the annual average of 2 m air temperature of North-Hemisphere over land as described in [75]. The GEV distribution is then estimated using the maximum likelihood method and we compute the upper bound using the formula (A.2).

We then report the 'upper' upper bound for $T$ in 2021:

$$z_{\max}(2021) = \mu_0 + \mu_1 T_{2021} - \frac{\sigma}{\xi}. \quad (A.3)$$

We compute the uncertainty range over this upper bound using $n=1000$ bootstrap samples of

the initial series. For each bootstrap series, we fit a non-stationary GEV law (with $T$ as a covariate on $\mu$) and compute the upper bound. We take the median values of these bootstrapped estimates as the best estimate and we define the lower (resp. higher) bound of the estimate as the 97.5% and 2.5% quantiles of the bootstrapped $z_{\max}$.

## Appendix B. Uncertainty in $T_{s,max}$ estimation

For computing the uncertainty range of the absolute $T_{s,max}$ for the urban areas, we use for $T_{500}$ and $z_{500}$ the spread of the 10 ensemble members of the ERA5 reanalysis simulation. The lowest (resp. highest) bound of the range uses the lowest (resp. highest) $T_{500}$ and $z_{500}$ over the members. The ensemble spread is small therefore the uncertainty range for the absolute $T_{s,max}$ is around 1 K.

For the computation of the more realistic $T_{s,max}$, we also need to take into account the uncertainty in the $q_{\min}$ value. To do so, we fit a stationary GEV law over the JJAS annual minimum of hourly 2 m air specific humidity over the 82 years between 1940–2021. This allows to give an approximate value to $P(q) := \mathbb{P}(q_{\min} \ge q)$ where $q$ is a fixed threshold. We compute the uncertainty range over this estimated probability using $n=1000$ bootstrap samples of the initial series. For each bootstrap series, we fit a GEV law and compute the return level $q_{82}$ associated with a probability $P(q_{82}) = \frac{1}{82}$, similar to the probability of the minimum observed in the ERA5 dataset. We then take the 97.5% and 2.5% quantiles of the bootstrapped $q_{82}$ as the uncertainty range for $q_{\min}$. Due to the resampling procedure, the 2.5% quantile is negative for all cities. As this is not physical, we take this lower bound to be null. We then combine the lowest (resp. highest) $T_{500}$ and $z_{500}$ over the members to the highest (resp. lowest) $q_{\min}$ to get the uncertainty range. Due to the large uncertainty in $q_{\min}$, the uncertainty range for the realistic $T_{s,max}$ is also large.
## Appendix C. Additional figures

<figure><img src="image_10.png" /></figure>

**Figure C1.** Analogs averaged daily mean 2 m air temperature T2M for (a) the counterfactual and (b) the factual periods. (c) Difference between the two periods of the averaged field. In panel (c), only the grid points with a significant difference at the 5% level are displayed.

<figure><img src="image_11.png" /></figure>

**Figure C2.** Vertical profiles of analogs averaged air temperature for (a) the Paris region (lat = [47.5°,49.5°], lon = [1°,5°]) and (c) the Netherlands region (lat = [51°,53.25°], lon = [5°,9°]). Vertical profiles of analogs averaged air specific humidity for (b) the Paris region and (d) the Netherlands region. The white and magenta dots show the mean lifting condensation level in the two periods. The associated vertical bars display a two standard deviation spread over the analogs.

<figure><img src="image_12.png" /></figure>

**Figure C3.** Analogs averaged daily mean upper layer soil moisture (10 cm) for (a) the counterfactual and (b) the factual periods. (c) Difference between the two periods of the averaged field. In panel (c), the hashed grid points correspond to a significant difference at the 5% level.
<figure><img src="image_13.png" /></figure>

**Figure C4.** June-July-August-September average of daily minimum of hourly mean specific humidity at 2 m for (a) the counterfactual period and (b) the factual period. (c) Difference between the two periods of the averaged field. In panel (c), the colored grid points represent a significant difference at the 5% level.

## References

[1] Pörtner H O *et al* 2022 IPCC sixth assessment report

[2] Seneviratne S I, Zhang X, Adnan M, Badi W, Dereczynski C, Di Luca A, Vicente-Serrano S M, Wehner M and Zhou B 2021 Weather and climate extreme events in a changing climate *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change* (Cambridge University Press)

[3] Campbell S, Remenyi T A, White C J and Johnston F H 2018 *Health Place* **53** 210–8

[4] Breshears D D, Fontaine J B, Ruthrof K X, Field J P, Feng X, Burger J R, Law D J, Kala J and Hardy G E S J 2021 *New Phytol.* **231** 32–39

[5] Huang W T K, Braithwaite I, Charlton-Perez A, Sarran C and Sun T 2022 *Environ. Res. Lett.* **17** 034017

[6] Masselot P *et al* 2023 *Lancet Planet. Health* **7** e271–81

[7] Hatfield J L and Prueger J H 2015 *Weather Clim. Extremes* **10** 4–10

[8] Brás T A, Seixas J, Carvalhais N and Jägermeyr J 2021 *Environ. Res. Lett.* **16** 065012

[9] Bastos A *et al* 2021 *Earth Syst. Dyn.* **12** 1015–35

[10] Zuo J, Pullen S, Palmer J, Bennetts H, Chileshe N and Ma T 2015 *J. Cleaner Prod.* **92** 1–12

[11] Barriopedro D, García-Herrera R, Ordóñez C, Miralles D and Salcedo-Sanz S 2023 *Rev. Geophys.* **61** e2022RG000780

[12] Perkins S E 2015 *Atmos. Res.* **164** 242–67

[13] Horton R M, Mankin J S, Lesk C, Coffel E and Raymond C 2016 *Curr. Clim. Change Rep.* **2** 242–59

[14] Domeisen D I, Eltahir E A, Fischer E M, Knutti R, Perkins-Kirkpatrick S E, Schär C, Seneviratne S I, Weisheimer A and Wernli H 2022 *Nat. Rev. Earth Environ.* **4** 1–15

[15] Petoukhov V, Rahmstorf S, Petri S and Schellnhuber H J 2013 *Proc. Natl Acad. Sci.* **110** 5336–41

[16] Coumou D, Petoukhov V, Rahmstorf S, Petri S and Schellnhuber H J 2014 *Proc. Natl Acad. Sci.* **111** 12331–6

[17] Kornhuber K, Coumou D, Vogel E, Lesk C, Donges J F, Lehmann J and Horton R M 2020 *Nat. Clim. Change* **10** 48–53

[18] Di Capua G, Sparrow S, Kornhuber K, Rousi E, Osprey S, Wallom D, van den Hurk B and Coumou D 2021 *npj Clim. Atmos. Sci.* **4** 1–14

[19] Seneviratne S I, Lüthi D, Litschi M and Schär C 2006 *Nature* **443** 205–9

[20] Fischer E M, Seneviratne S I, Lüthi D and Schär C 2007 *Geophys. Res. Lett.* **34** L05707

[21] Seneviratne S I, Corti T, Davin E L, Hirschi M, Jaeger E B, Lehner I, Orlowsky B and Teuling A J 2010 *Earth-Sci. Rev.* **99** 125–61

[22] Hirschi M, Seneviratne S I, Alexandrov V, Boberg F, Boroneant C, Christensen O B, Formayer H, Orlowsky B and Stepanek P 2011 *Nat. Geosci.* **4** 17–21

[23] Miralles D G, Van Den Berg M, Teuling A and De Jeu R 2012 *Geophys. Res. Lett.* **39** L21707

[24] Miralles D G, Teuling A J, Van Heerwaarden C C and Vilà-Guerau de Arellano J 2014 *Nat. Geosci.* **7** 345–9

[25] Rasmijn L, Van der Schrier G, Bintanja R, Barkmeijer J, Sterl A and Hazeleger W 2018 *Nat. Clim. Change* **8** 381–5

[26] Dirmeyer P A, Balsamo G, Blyth E M, Morrison R and Cooper H M 2021 *AGU Adv.* **2** e2020AV000283

[27] Petoukhov V, Petri S, Rahmstorf S, Coumou D, Kornhuber K and Schellnhuber H J 2016 *Proc. Natl Acad. Sci.* **113** 6862–7

[28] Kornhuber K, Petoukhov V, Karoly D, Petri S, Rahmstorf S and Coumou D 2017 *J. Clim.* **30** 6133–50

[29] Mann M E, Rahmstorf S, Kornhuber K, Steinman B A, Miller S K and Coumou D 2017 *Sci. Rep.* **7** 1–12

[30] Mann M E, Rahmstorf S, Kornhuber K, Steinman B A, Miller S K, Petri S and Coumou D 2018 *Sci. Adv.* **4** eaat3272

[31] Kornhuber K and Tamarin-Brodsky T 2021 *Geophys. Res. Lett.* **48** e2020GL091603

[32] Ghil M *et al* 2011 *Nonlinear Process. Geophys.* **18** 295–350

[33] Fisher R A and Tippett L H C 1928 Limiting forms of the frequency distribution of the largest or smallest member of a sample *Mathematical Proceedings of the Cambridge Philosophical Society* vol 24 (Cambridge University Press) pp 180–90

[34] Gnedenko B 1943 *Ann. Math.* **44** 423–53

[35] Coles S, Bawa J, Trenner L and Dorazio P 2001 *An Introduction to Statistical Modeling of Extreme Values* vol 208 (Springer)

[36] Cheng L, AghaKouchak A, Gilleland E and Katz R W 2014 *Clim. Change* **127** 353–69

[37] National Academies of Sciences Engineering and Medicine 2016 *Attribution of Extreme Weather Events in the Context of Climate Change* (The National Academies Press)

[38] Philip S *et al* 2020 A protocol for probabilistic extreme event attribution analyses *Adv. Stat. Climatol. Meteorol. Oceanogr.* **6** 177–203

[39] Van Oldenborgh G J *et al* 2021 *Clim. Change* **166** 13

[40] Van Oldenborgh G J, Wehner M F, Vautard R, Otto F E, Seneviratne S I, Stott P A, Hegerl G C, Philip S Y and Kew S F 2022 *Earth's Future* **10** e2021EF002271

[41] Auld G, Hegerl G and Papastathopoulos I 2021 Adv. Stat. Clim. Meteorol. Oceanogr. **9** 45–66

[42] Mitchell D, Kornhuber K, Huntingford C and Uhe P 2019 *Lancet Planet. Health* **3** e290–2

[43] Philip S Y *et al* 2022 *Earth Syst. Dyn.* **13** 1689–713

[44] Thompson V, Kennedy-Asser A T, Vosper E, Lo Y E, Huntingford C, Andrews O, Collins M, Hegerl G C and Mitchell D 2022 *Sci. Adv.* **8** eabm6860

[45] Vautard R *et al* 2020 *Environ. Res. Lett.* **15** 094077
[46] Patterson M 2023 *Geophys. Res. Lett.* 50 e2023GL102757

[47] Wouters J and Bouchet F 2016 *J. Phys. A: Math. Theor.* 49 374002

[48] Ragone F, Wouters J and Bouchet F 2018 *Proc. Natl Acad. Sci.* 115 24–29

[49] Yiou P and Jézéquel A 2020 *Geosci. Model Dev.* 13 763–81

[50] Ragone F and Bouchet F 2021 *Geophys. Res. Lett.* 48 e2020GL091197

[51] Gessner C, Fischer E M, Beyerle U and Knutti R 2021 *J. Clim.* 34 6619–34

[52] Zhang Y and Boos W R 2023 *Proc. Natl Acad. Sci.* 120 e2215278120

[53] Yiou P, Jézéquel A, Naveau P, Otto F E, Vautard R and Vrac M 2017 *Adv. Stat. Climatol. Meteorol. Oceanogr.* 3 17–31

[54] Jézéquel A, Yiou P and Radanovics S 2018 *Clim. Dyn.* 50 1145–59

[55] Faranda D, Bourdin S, Ginesta M, Krouma M, Messori G, Noyelle R, Pons F and Yiou P 2022 *Weather Clim. Dyn.* 3 1311–40

[56] Sousa P M, Barriopedro D, García-Herrera R, Ordóñez C, Soares P M and Trigo R M 2020 *Commun. Earth Environ.* 1 48

[57] Ma F, Yuan X, Jiao Y and Ji P 2020 *Geophys. Res. Lett.* 47 e2020GL087809

[58] Hersbach H et al 2020 *Q. J. R. Meteorol. Soc.* 146 1999–2049

[59] Cadiou C, Noyelle R, Malhomme N and Faranda D 2023 *Asia Pac. J. Atmos. Sci.* 59 83–94

[60] Faranda D, Pascale S and Bulut B 2023 *Environ. Res. Lett.* 18 034030

[61] Xu P, Wang L, Huang P and Chen W 2021 *Atmos. Res.* 252 105446

[62] Maraun D 2016 *Curr. Clim. Change Rep.* 2 211–20

[63] Zhang L 2016 *Clim. Dyn.* 47 3157–69

[64] van Haren R, van Oldenborgh G J, Lenderink G and Hazeleger W 2013 *Environ. Res. Lett.* 8 014053

[65] Welch B L 1947 *Biometrika* 34 28–35

[66] Chadwick R, Good P and Willett K 2016 *J. Clim.* 29 7613–32

[67] Coffel E D, Horton R M, Winter J M and Mankin J S 2019 *Environ. Res. Lett.* 14 084003

[68] Padrón R S, Gudmundsson L, Decharme B, Ducharne A, Lawrence D M, Mao J, Peano D, Krinner G, Kim H and Seneviratne S I 2020 *Nat. Geosci.* 13 477–81

[69] McKinnon K A, Poppick A and Simpson I R 2021 *Nat. Clim. Change* 11 598–604

[70] Merrifield A L, Simpson I R, McKinnon K A, Sippel S, Xie S P and Deser C 2019 *Geophys. Res. Lett.* 46 14082–92

[71] Li M, Wu P and Ma Z 2020 *Int. J. Climatol.* 40 5744–66

[72] Schumacher D L, Keune J, Van Heerwaarden C C, Vilà-Guerau de Arellano J, Teuling A J and Miralles D G 2019 *Nat. Geosci.* 12 712–7

[73] Bieli M, Pfahl S and Wernli H 2015 *Q. J. R. Meteorol. Soc.* 141 98–108

[74] Sánchez-Benítez A, Goessling H, Pithan F, Semmler T and Jung T 2022 *J. Clim.* 35 2373–90

[75] Osborn T J, Jones P D, Lister D H, Morice C P, Simpson I R, Winn J, Hogan E and Harris I C 2021 *J. Geophys. Res.: Atmos.* 126 e2019JD032352