Distinct influences of large-scale circulation and regional feedbacks in two exceptional 2019 European heatwaves

Two separate heatwaves affected western Europe in June and July 2019, in particular France, Belgium, the Netherlands, western Germany and northeastern Spain. Here we compare the European 2019 summer temperatures to multi-proxy reconstructions of temperatures since 1500, and analyze the relative influence of synoptic conditions and soil-atmosphere feedbacks on both heatwave events. We find that a subtropical ridge was a common synoptic set-up to both heatwaves. However, whereas the June heatwave was mostly associated with warm advection of a Saharan air mass intrusion, land surface processes were relevant for the magnitude of the July heatwave. Enhanced radiative fluxes and precipitation reduction during early July added to the soil moisture deficit that had been initiated by the June heatwave. We show this deficit was larger than it would have been in the past decades, pointing to climate change imprint. We conclude that land-atmosphere feedbacks as well as remote influences through northward propagation of dryness contributed to the exceptional intensity of the July heatwave. Two exceptional heatwaves affected Europe in June and July 2019. Although both were driven by the large-scale circulation, the July event was also amplified by soil moisture feedback, according to an analysis of past temperatures, weather analogues and soil-atmosphere interactions.

H eatwaves (HWs) are among the most concerning extreme meteorological events, as they have a wide range of impacts, including human health (e.g. increased mortality and morbidity) 1,2 and significant socio-economic and ecological effects, such as wildfires and poor air quality events 3,4 , droughts 5 and peaks in energy consumption demand 6,7 . Within the context of global warming, an increased frequency in extremely warm events is foreseen, comprising HWs of unprecedented extension and duration [8][9][10] . 2019 was the second warmest year at the global scale, only surpassed by the strong El-Niño year of 2016 11 . Unsurprisingly, summer 2019 presented exceptional HWs in Europe, exceeding notorious episodes which occurred just 1 year before in the also very hot summer of 2018 12,13 . In terms of affected areas, the 2019 HW events resembled to a large extent the 2003 summer HW 14,15 and in many places temperature extremes even shattered those of 2003. In late June an outstanding HW began in southwestern Europe 12 , and extended towards most of France and parts of central Europe. During this event the city of Vérargues in southeastern France reached an astonishing daily maximum temperature (TX hereon) of 46°C on June 28th. This was the first time temperature measurements exceeded 45°C in France. Just a few weeks later, another exceptional HW set new historical values in France and other European countries. For example, Paris registered a TX of 43°C, surpassing the previous record standing since 1947 by~2°C. Furthermore, for the first time since the beginning of meteorological observations, Belgium and the Netherlands exceeded the 40°C barrier. Fortunately, summer 2019 caused considerably less mortality excess than previous HWs, including the devastating 2003 event 16 . This might result from the combination of human factors, which include the lessons learned from the 2003 HW (i.e. early warning systems, better preparedness and societal awareness, deployment of sheltering and water-cooling facilities, use of air conditioning, etc.), and the shorter duration of both 2019 HWs.
Most extreme temperature events are partially driven by anomalous large-scale atmospheric circulation. However, the current rate of warming (i.e. thermodynamic changes) is sufficient to produce exceptional HWs, even without unprecedented anomalies in the large-scale circulation. Contrasting with the lack of robust projections in dynamical changes 17 , recent works indicate robust and significant increases in maximum HW magnitude over large regions, even for 1.5°C global warming targets 8 . Moreover, anthropogenic forcing has already caused a 7fold increase in the likelihood of extreme heat events 18 . In addition to direct radiative effects of increasing greenhouse gases concentrations, the potential contribution of enhanced local land-atmosphere feedbacks has also been acknowledged 19 . Recent studies have further explored non-local feedbacks in recent mega-HW events 20,21 . Local drying and subsequent enhanced surface heat fluxes, together with horizontal warm advection and heat accumulation in the atmospheric boundary layer, have been shown to contribute to the magnitude of temperature anomalies during the August 2003 HW 22 . Very similar processes were observed in the 2010 Russian mega-HW 22,23 . Recently, the same events have been explored to introduce the concept of upwind soil dryness 21 (also referred to as selfpropagation 20 ). These conceptual models illustrate how air masses warmed by sensible heat fluxes (due to pre-conditioning dry soil conditions) can be advected downwind, stimulating land-atmosphere feedbacks in nearby regions that contribute to a progressive set-up for HW occurrence.
Nevertheless, the combined effect of regional scale feedback processes and large-scale atmospheric circulation is crucial for understanding the development of extreme heat events 24 . Regarding the latter, several studies have described the weather systems associated with European HWs, highlighting the key role of blocking/ridge occurrence 12,[25][26][27] . Given this major control of the dynamics, other studies have used flow-analogue approaches to quantify the contribution of thermodynamic changes to the observed magnitude of outstanding recent HW events [28][29][30] . Here we aim to investigate the combined roles of large-scale atmospheric circulation and soil moisture conditions in the occurrence and severity of the summer 2019 HWs in Europe. The synoptic conditions are analyzed for both events to characterize large-scale circulation features. Thermodynamic aspects are also considered, including regional feedback processes and their lagged effects. In particular, we investigate soil-atmosphere processes observed since the June 2019 HW and their potential role in amplifying the magnitude of the July 2019 HW.

Results
How anomalous was summer 2019 in Europe? To place the 2019 European summer into a long-term context, we first derived estimates of the European mean summer temperature anomalies since 1500 using a multi-proxy reconstruction (1500-1900) 31 and an observational dataset (1901-2019) 32 . Further details are provided in the Data and Methods section. Summer 2019 falls within the top five warmest recorded in Europe since the early 16th century, only surpassed by other recent devastating summers (2018, 2010 and 2003; Fig. 1a). Although extremely warm conditions for the 2019 summer were mainly confined to western and some central areas of the continent, temperature anomalies were large enough to produce an overall continental anomaly close to +2°C (with respect to 1981-2010). Recent results estimate a return period of nearly 300 years for a similar event, taking into account recent climatic conditions 18 . As seen in Fig. 1b, warm summers have become more frequent since the last decades of the 20th century, and the 21st century concentrates an unusual frequency of extreme summers when compared to the long-term variability (1500-2019). This is reflected by the dominance of post-2000 events in the high-end tail of the European summer temperature distribution (histogram, Fig. 1a), and the pronounced shift in the 30-year Gaussian fitted distributions between 1960-1989 and 1990-2019 (Fig. 1a, light and dark grey shading).
Interestingly, the European mean temperature anomaly of summer 2019 falls very close to that of summer of 2003 (Fig. 1a). The spatial distribution of TX anomalies averaged for summer 2019 was also somewhat similar to that observed in 2003, according to the E-OBS dataset (see Fig. 1c). Taking into account the comparable magnitude and spatial signatures of the 2003 14,33 and 2019 summers, we used the former as a "benchmark" to evaluate the exceptionality of the 2019 summer temperature anomalies. Given the fast pace of current atmospheric warming, the comparison of 2003 and 2019 summers was performed by defining temperature anomalies with respect to two distinct baselines: (i) the full available period (1950-2018) and (ii) the previous 30-year period at the time these summers occurred (i.e. 1973-2002 for 2003 and 1989-2018 for 2019), which leads to a warmer climatological baseline for 2019 than for 2003. As shown in Fig. 1c, summer mean TX anomalies were in general larger for 2003 than in 2019 with respect to their corresponding climatological conditions. However, the highest daily TX anomalies registered in 2019 exceeded those observed in 2003, even if anomalies are computed with respect to their corresponding climatologies (Fig. 1d). Specifically, daily TX in northeastern France and Benelux surpassed climatological values by nearly 20°C during summer 2019, compared with the maximum TX anomalies observed in the same regions during summer 2003 (~16°C). Note that these values are also the warmest TX anomalies of the continent in both summers. Therefore, the dichotomy between summer averages and daily values reflects that while summer 2003 was more anomalous for the climatic conditions expected at that time (essentially due to the longer nature of the August 2003 HW), the 2019 HWs were generally more intense than the 2003 HW at daily time scales in many places of western and central Europe. Similar conclusions are obtained when using the full period (1950-2018) as baseline ( Supplementary Fig. 1). To illustrate the contribution of the non-stationary climate conditions to the magnitude of recent recordbreaking events, Supplementary Fig. 2a shows the difference between the summer mean TX 30-year climatologies as of 2019 and 2003 in Europe. In just a~15-year period, TX normals have increased by more than 1°C in summer over most areas (and even by~2°C at some locations of southern Europe). Additionally, the rate of record-breaking events has also been increasing To stress the exceptionality of the 2019 HWs, Fig. 2a presents the areas where all-time records in TX were hit during that summer. The spatial pattern shows a strong resemblance with the corresponding map for summer 2003, in particular over France and Benelux (see Supplementary Fig. 3). This means that a substantial part of the records set in summer 2003 was broken during 2019, with the exception of some areas in southwestern Europe (e.g. western Iberia, where 2003 temperatures were shattered during summer 2018 12 ).
The unprecedented temperatures during summer 2019 were associated with two clearly distinct HWs in late June and late July. In Fig. 2b, c the spatial distribution and duration of these HW events are depicted, as diagnosed from a novel HW tracking algorithm (see "Data and Methods"). The panels show that the spatial distribution of areas under HW conditions during July extended much further north than those during the June HW. This is in agreement with the timing of new all-time records presented in Fig. 2a. During July, unprecedented TX was reported over larger areas and dominated higher latitudes, including more than half of the French territory, the Benelux, western Germany, southeastern England and parts of Scandinavia. In contrast, daily all-time records in June were essentially restricted to southeastern France and northeastern Iberia. In spite of this, the highest absolute values (TX > 45°C) were observed during the June HW. The persistence of HW conditions was also more prominent over land during the June HW (cf. dots in Fig. 2b, c), with large areas of western Europe experiencing an extremely high number of HW days. These differences in the spatial signatures of the 2019 HWs are also reflected in the latitudinal location of the 500 hPa geopotential height (Z500) anomalies (contour lines in Fig. 2b, c), suggesting distinct atmospheric circulation patterns.
Atmospheric circulation during the 2019 HWs. In this section, we describe the large-scale atmospheric circulation configurations behind the summer 2019 HWs. The temporal evolution and spatial tracking of the two HWs (summarized in Supplementary  Fig. 4) show that the initial location of the HW centre was detected much further south in June than in July. In the former, HW conditions originated over northern Africa and then migrated towards northern France, before affecting eastern Europe during  its later stages. In contrast, the July HW onset was detected over France and then moved to higher latitudes, reaching areas close to the Arctic towards the end of its lifecycle.
Despite these differences, a relevant common factor can be identified. Both events displayed a classical pattern of Z500 positive anomalies over the affected areas, accompanied by the presence of a low-pressure system in the eastern Atlantic (see also Supplementary Fig. 4). 1000-500 hPa geopotential height thicknesses averaged for the HW periods are presented in Fig. 3a, b, revealing pronounced ridge-like patterns in both events, extending from northern Africa towards western Europe. However, the ridge affecting southwestern Europe during the June HW was stronger and better defined (i.e. sharper zonal gradients). As a result, a stronger southerly wind component Fig. 3 Synoptic configuration and forcing mechanisms. a Number of days (shading) when a Saharan intrusion was detected in each grid cell during the June HW. Black dots represent areas where the occurrence of Saharan intrusions was unprecedented (since at least 1948). Lines depict the composite of 1000-500 hPa geopotential height thickness (dam) for the days when our algorithm detects heatwave conditions. The dashed contour at 580 dam indicates the minimum threshold for Saharan dust intrusion. Panel adapted from Sousa et al. 12 . b Same as (a), but for the July HW. c Temporal evolution of the 1000-850 hPa temperature anomaly (°C, with respect to 1981-2010; grey line; right y-axis) and the contributors to the temperature tendency (°C day −1 ; coloured lines; left y-axis) over WEU, from lag −8 to lag +8 days of the June HW onset there (24 June). Grey shading represents days with HW conditions in WEU. d Same as (c) but for the July HW (onset on 23 July). e Temporal evolution of the fractional area (%) dominated by each forcing during the June HW. Grey shading as for previous panels. f Same as (e) but for the July HW. A 3-day smoothing is applied to the series presented in (c-f).
COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-020-00048-9 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | (2020) 1:48 | https://doi.org/10.1038/s43247-020-00048-9 | www.nature.com/commsenv characterized this first HW, when compared to relatively more stagnant conditions during the late July episode. This is in agreement with the strong intensity of the Saharan intrusion (see "Data and Methods" for details) observed during the first event ( Fig. 3a), when an air mass with desertic features reached unprecedented latitudes over France. Saharan intrusions have been shown to be associated with extreme heat events in southwestern Europe 12,34,35 , as they present very high potential temperatures and low moisture content, favoring intense surface warming under anticyclonic conditions. During the July HW, air masses with such thermodynamic properties were not detected further north than the western Mediterranean (Fig. 3b). In consequence, these results (supported by the HWs evolution) suggest a more pronounced influence of warm advection during the June HW.
Following the Eulerian methodology developed in previous studies 12,26 , in Fig. 3c, d the main physical mechanisms contributing to the temperature anomalies in the lower troposphere are examined for both HWs (see "Data and Methods"). Air motions, both horizontal (warm advection, red line) and vertical (strong adiabatic heating due to subsidence, blue line) are often important for the establishment and maintenance of HWs over Europe, although their relative contributions may differ 22,36 . This is also the case for the 2019 HWs. For the June HW, horizontal advection (red) was relevant before and at the onset of the HW, while vertical descent (blue) was essential for its maintenance (Fig. 3c). Differently, diabatic processes (green line) played a more important role for the setup of the July HW (Fig. 3d). This diversity in the underlying processes of the HW events is even more noticeable considering the fraction of areas (within western Europe, WEU, [43°-53°N, 0°-10°E]) where each contributing factor accounted for the largest temperature changes during the HWs lifecycles (Fig. 3e, f). Accordingly, as discussed in detail below, regional diabatic processes played a key role during the July HW. In Supplementary Fig. 5 (top panels) these differences are reinforced by the day-to-day evolution of the vertical profiles of temperature anomalies and horizontal wind averaged over WEU. During the June HW air masses presented reduced vertical gradients (presumably associated with the presence of the vertically homogeneous Saharan warm air intrusion) as compared to the July HW. Also, wind vectors support the major role of horizontal advection in the onset of the June HW, contrasting with more stagnant conditions at the peak of both HWs. This is further evidenced by the mean WEU vertical profiles of absolute temperature averaged over the HWs duration, as well as the instantaneous profiles at the peak of the HWs (bottom panels of Supplementary Fig. 5). Moreover, during the July HW, temperature anomalies seem to propagate upwards from the surface a few days prior to the HW onset, suggesting a progressive surface-atmosphere coupling, building up in the lower troposphere towards the HW peak. A similarly gradual warming process has been reported in the boundary layer, leading to self-intensification of near-surface temperatures during the well-known mega-HWs observed in western (2003) and eastern (2010) Europe 22 . In the next section, we further explore whether the diabatic processes that dominated the establishment of the July HW were influenced by land-atmosphere feedbacks.
Amplification of the late July 2019 HW due to soil desiccation. To explore the presence of land-atmosphere feedbacks during the July HW, we first analyzed the temporal evolution of a set of relevant variables averaged over a box covering the region with the highest TX anomalies (northeastern France and Belgium) as presented in Fig. 4. The preceding June HW contributed to strong losses in soil moisture content in that region, with this drying also reinforced by above-average radiative fluxes at the surface and low precipitation throughout July (see Supplementary Fig. 6). Consequently, persistent soil desiccation occurred between the two HWs. This resulted in anomalous surface heat fluxes, as shown in the lower panel of Fig. 4a. The three-week period in between the two HWs was characterized by an approximate doubling (halving) of the sensible (latent) heat fluxes when compared to the corresponding climatological values for that time of the year. This is reflected in the recurrence of days with Bowen ratio values above 1, indicating that energy partition was dominated by sensible heat fluxes from the surface, due to soil moisture limited latent fluxes (see also Supplementary Fig. 6). These results point to a contribution from regional soil moisture deficit to near-surface warming that persisted until the onset of the July HW (i.e. local land-atmosphere processes). Figure 4b illustrates relevant fields for the land-atmosphere coupling on a larger spatial domain than the regional box over NE France/Belgium considered in Fig. 4a. The areas with significant soil dryness (dots) in the weeks preceding the July HW are in good spatial agreement with subsequent large sensible heat flux positive anomalies (shading) during the build-up of the July HW over NE France/Belgium. Collocated large anomalies of both fields extended over large areas, suggesting land-atmosphere coupling beyond NE France/Belgium, particularly to the south of the region hit by the July HW. This, along with the mean nearsurface wind direction observed in the days preceding the July HW, suggests similar processes to those describing the concept of self-propagation 21,22 . Under the presence of southerly winds, dry air masses in central France warmed by anomalous sensible heat fluxes prior to the July HW were advected further north, likely enhancing remote land-atmosphere feedbacks and local sensible heat fluxes over the box displayed in Fig. 4b. This process arguably contributed to the amplification of the July HW. The circulation analogue exercise conducted further ahead supports these conclusions.
A similar analysis was performed for the region with the highest temperature anomalies during the June HW (southeastern France, Supplementary Fig. 7). During the intense June HW, energy transfer from the surface by sensible heat fluxes was below climatological values. In addition, in this area close to the Mediterranean, soil desiccation was not as intense as observed further north. These results suggest land-atmosphere coupling did not substantially contribute to the June HW, thus reinforcing the contrast between the two events, i.e. the more advective nature of the June HW compared to the dominance of diabatic processes associated with land-atmosphere coupling during the July HW.
To further deepen the process analysis discussed above, Figs. 5 and 6 present two distinct analogue exercises (see "Data and Methods") with the aim of evaluating: (i) the potential contribution of the June HW to the subsequent soil desiccation observed in July; (ii) the level of amplification of surface temperature anomalies during the July HW as a result of the preceding soil moisture deficits.
The results of the flow analogues for the June HW indicate that recent circulation conditions similar to those reported during the June HW have some drying imprints in the subsequent soil moisture conditions of western Europe (Fig. 5b). Indeed, the soil moisture content over WEU (dashed box in Fig. 5c) is significantly lower for flow analogues of the June HW (Fig. 5d, dark boxes) than for random circulation conditions (light boxes; p < 0.05; t-test and Kolmogorov-Smirnov test), portraying the role of the atmospheric circulation pattern in driving subsequent soil moisture deficits. These differences are even larger when using ERA5 (1979-2019) or ERA20C (1900-2010) data as a pool of analogues (see Supplementary Figs. 8 and 10), suggesting reduced variability of soil moisture in NCEP/NCAR (i.e. weaker responses to atmospheric forcing) and/or differences related to the thickness of the uppermost soil layer (0-10 cm in NCEP/ NCAR vs. 0-7 cm in ERA reanalyses). These results lend support to the hypothesis that the atmospheric circulation associated with the June HW contributed to the subsequent desiccation that preceded the July HW. Note that drying was not so obvious during the actual June HW in observations (Fig. 4a), arguably because soils were replenished throughout a deeper layer (0-2 m) prior to the event (upper panel of Supplementary Fig. 6), which might have contributed to an initial dampening of the desiccation process. On the other hand, the results of the analogue exercise also indicate that flow analogues precede drier conditions in the present than in the recent past (Fig. 5a, b). Part of this difference is associated with a generalized regional drying over the analyzed period, since a comparable soil moisture decrease is also observed between the random circulation distributions of both subperiods (Fig. 5d), which are not constrained by the atmospheric circulation. Qualitatively similar results are obtained for ERA reanalyses (see Supplementary Figs. 8d and 10a), although trends and patterns are overall weaker in ERA20C, arguably due to the lack of soil moisture-related observations in the assimilation process. The temporal differences between soil moisture distributions are consistent with the reported occurrence of more severe European droughts due to enhanced atmospheric evaporative demand by recent warming trends 37,38 . In summary, our results support that the June HW, together with the precipitation deficits and high radiative fluxes that followed it, contributed to the soil moisture deficits preceding the July HW. Furthermore, we also find that this drying signal has been amplified in recent decades. While this result should not be interpreted as a formal attribution to anthropogenic factors, it is in agreement with recent studies attributing dry-season water imbalance changes to humaninduced climate change 39 .
To support the above-mentioned amplifying role of the observed soil moisture deficits in the magnitude of the July HW, we have searched for flow analogues of each day of the July HW and reconstructed the associated TX anomalies (Fig. 6). We account for the role of soil desiccation by distinguishing between analogue days preceded by dry and wet conditions over WEU, as inferred from regional mean anomalies of soil moisture averaged for the previous 15 days (see "Data and Methods"). The results indicate that similar flow patterns to those recorded during the July HW tend to cause warmer conditions when they are preceded by dry conditions (Fig. 6a, b). In other words, for similar atmospheric circulation, soil moisture deficits promote warming (Fig. 6d, dark boxes; see also Fig. 6c). By construction, this warming should be interpreted as a response to drying, and not averaged during the 7 days prior to the July HW, with vectors depicting the mean near-surface wind during the same period. Dots represent areas with soil moisture deficits (different sizes depict 10%, 20 and 30% deficit), averaged during the 15 days before the July HW onset (here using the C3S satellite derived product). Black box represents the area considered for the series presented in a).
COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-020-00048-9 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | (2020) 1:48 | https://doi.org/10.1038/s43247-020-00048-9 | www.nature.com/commsenv the other way around, since trends have been removed and the use of time lags minimizes misattributions of cause and effect. Random distributions (unconstrained by the atmospheric circulation) indicate similar warming levels following short-term soil moisture deficits (Fig. 6d, light boxes). Accordingly, regional drying seems to favour above-normal temperatures, regardless of the atmospheric circulation. Interestingly, additional analyses reveal atmospheric circulation differences between the flow analogues of dry and wet years (Fig. 6e, dark whiskers), involving larger positive Z500 anomalies for flow analogues preceded by soil moisture deficits (see Fig. 6c), which translate to lower RMSE (i.e. closer patterns to the actual circulation) than during wet conditions. These differences in RMSE are also observed for the unconditional distributions (Fig. 6e, light whiskers), indicating an overall tendency for dry periods to precede higher pressure anomalies. This could reflect methodological issues (e.g. limited sampling, residual trends, autocorrelation issues), although we obtain similar results for dry and wet periods of the ERA5 (1979-2019) and ERA20C (1900-2010) reanalyses (see Supplementary Figs. 9 and 10). Alternatively, the Z500 differences between dry and wet conditions may also indicate feedbacks of soil moisture deficits on the atmospheric circulation anomalies. If the latter is the case, such effect herein involves somehow weak high-pressure patterns, whose spatial details depend on the considered dataset (cf. Fig. 6c and Supplementary Fig. 9c) and methodological choices. Previous studies have suggested atmospheric circulation responses to soil moisture deficits, including local effects through thermal expansion by enhanced sensible heat fluxes 40 , and remote effects caused by a thermally-induced low 41 or changes in cloud cover 42 . In short, our results indicate that soil moisture deficits in western Europe intensified the warming already expected from the circulation observed during the July HW and might even have contributed to amplifying the circulation anomalies. Additional studies are warranted to explore and quantify the contribution of atmospheric circulation responses induced by land-atmosphere coupling to the intensity of HWs.
Summary and discussion. Two distinct HWs affected widespread areas of western Europe in June and July 2019, contributing to placing that summer within the top five warmest since 1500 at the European scale. While the spatial distribution of the affected areas strongly resembled the historical HW of August 2003, the relatively shorter-lived 2019 HWs were more intense on daily time scales, shattering previous all-time records in many places (some of them standing since 2003). Here we have dissected these events with recently developed tools to provide an assessment of different relevant factors: (i) the role of the dynamics (synoptic setups associated with the 2019 HWs); (ii) underlying physical processes (warm advection vs. diabatic fluxes, including enhanced near-surface heating due to soil moisture deficits); (iii) recent thermodynamic changes (the steady regional warming trend).
We have applied recent novel methodologies to track the two 2019 HWs and the occurrence of subtropical warm air intrusions. The June HW displayed a clear fingerprint of a Saharan intrusion. The advection of exceptionally warm and dry air, together with enhanced subsidence under pronounced Z500 anomalies, are the main features of this event that triggered all-time temperature records in southern France and northeastern Spain. Differently, the July HW extended much further north, and unprecedented temperatures hit a comparatively larger domain of Europe. Diabatic processes, rather than temperature advection associated with a Saharan intrusion, played a dominant role for the setup of this event. Some studies have found different relative contributions of horizontal advection, subsidence and diabatic processes in shaping European HWs 22,36 . Our analysis attests to the distinct nature of two HWs that took place over the same region within a few weeks, supporting the coexistence of distinct dominant forcing mechanisms for their onset and maintenance.
We have shown evidence supporting a contribution of landatmosphere coupling to the temperature anomalies during the July HW, potentially involving the self-propagation mechanism discussed for previous HWs 20, 21 . The atmospheric conditions prevailing since the onset of the June HW, i.e. prolonged periods with no rain and persistently high solar radiative fluxes and temperatures, significantly contributed to anomalous soil moisture deficits over large areas, in particular over France during the transition period between both HWs. This resulted in strong energy transfer between the soil and atmosphere, via increased (decreased) sensible (latent) heat fluxes before the onset of the July HW. These diagnostics of land-atmosphere feedbacks were not restricted to the region most affected by the July HW (northeastern France and Belgium). They were also observed further south in areas under the influence of sustained southerlies, which suggests a northward propagation of dryness through the advection of warm air masses. By using atmospheric flow analogues, our analysis further supports the role of soil desiccation on the amplification of the July 2019 HW. In this context, lessened soil-atmosphere feedbacks have been reported in areas where shallow groundwater is available 43 . Accordingly, we argue that an event occurring under similar synoptic patterns to those observed in July 2019 would result in lower temperature anomalies if preceding soil conditions were wetter. Our results also suggest non-negligible effects of the June HW on the July HW through an imprint of the former on the soil moisture deficits that influenced the latter. However, further modelling studies are warranted to support this conclusion, as well as to address land-atmosphere feedbacks on the atmospheric circulation. In this regard, recent ensemble experiments provide evidence about complex local and remote effects of soil moisture deficits up to two months, including non-local responses in atmospheric circulation that can further amplify HW magnitudes 42 .
Our results also indicate that previous colder climatic conditions would have resulted in less soil desiccation than observed. This effect is found regardless of the atmospheric circulation, pointing to atmospheric warming effects as a consequence of anthropogenic forcing 44 , probably enhanced by land-use and land-cover changes 43,45 . We acknowledge limitations (e.g. limited sample size, transient climate conditions over the analyzed period or biases in the reanalysis datasets) and assumptions (e.g. event definition) in our analogue experiments. Moreover, their results should not be interpreted as a formal attribution to anthropogenic forcing, despite the overall consistency based on independent evidence. For example, following earlier studies on European HWs 46 , recent findings also indicate that hot summer temperatures in the Fig. 6 Soil-atmosphere feedbacks during the July HW. Mean Z500 (m, contours) and TX (°C, shading) anomalies (with respect to 1981-2010) reconstructed for the July HW (23-26 July 2019) from daily flow analogues of Z500 over Europe (solid box in c) preceded by (a) wet (above 66th) and (b) dry (below 33rd) soil moisture conditions at 0-10 cm in WEU (dashed box in c) during the previous 15 days. c Difference between (b) and (a). d Flowconditioned (dark grey boxes) and random (light grey boxes) distributions of the mean TX anomalies for the July HW over WEU preceded by wet and dry conditions (x-axis). Horizontal dashed line depicts regional mean averaged values observed during the July 2019 HW. e As (d) but for the root-meansquare error (RMSE) distributions of the flow-conditioned (dark grey boxes) and random (light grey boxes) analogues. Boxes and whiskers show the ±0.5·SD around the mean and 5th-95th percentile ranges, respectively, with circles denoting the maximum and minimum values. Data source: NCEP/ NCAR reanalysis.
Mediterranean area are often shortly preceded by the occurrence of dryness in spring or even early summer 47 , thus favoring temporally compounding events 48 . These facts highlight that future extreme heat episodes will probably be even further exacerbated by the increasing severity of drought events 37,49 , as stronger losses by evaporation are expected in a warming climate 50 whenever soil moisture is still available 51 .
Data and methods E-OBS dataset. Daily minimum and maximum 2 m temperature from the E-OBS gridded dataset (v21.0) was used to characterize anomalies and extremes during the 2003 and 2019 events, as well as trends and record-breaking values during the available period (since 1950). Anomalies are computed by removing the daily climatological mean (1981-2010). E-OBS is a European land-only high-resolution gridded observational dataset, using the European Climate Assessment and Dataset (ECA&D) blended daily station data 52 . It is presented on a horizontal resolution of 0.25°×0.25°. Files are replaced in monthly updates and in updated versions of the E-OBS dataset. Accordingly, small changes might occur between these releases after new data and/or stations are added.
NCEP/NCAR dataset. Meteorological fields were retrieved from the NCEP/NCAR reanalysis daily dataset 53 , starting from 1948. The following variables were considered for pressure levels between 1000 and 500 hPa on a 2.5°× 2.5°horizontal resolution grid: air temperature, geopotential height, zonal/meridional wind components, vertical velocity. We also analyzed other fields represented in a Gaussian grid: surface net radiation fluxes (longwave and short-wave), latent and sensible heat fluxes, precipitation, 2 m temperature, 10 m wind, potential evapotranspiration and soil moisture fraction (0-10 cm and 10-200 cm). These fields were used to: (i) characterize and track the HW events, (ii) derive a catalogue of Saharan intrusions, (iii) generate vertical profiles, (iv) compute the contributing terms to the temperature tendency equation, and (v) perform the analogue exercises. Specific methods for products derived from these variables are explained below. In all cases, anomalies are computed with respect to the climatological seasonal cycle (1981-2010).
ERA5 and ERA20C datasets. Meteorological fields were extracted from two ECMWF (European Centre of Medium-range Weather Forecast) reanalyses to replicate the analogue exercises (see methodology further ahead) performed with the NCEP/ NCAR dataset. The ERA5 54 and ERA20C 55 datasets were considered, using the highest horizontal resolution available for the latter (1.25°× 1.25°), for the 1979-2019 and 1900-2019 periods, respectively. Daily time series of 2 m temperature, Z500 and soil moisture fraction (0-7 cm) were retrieved.
European temperature reconstruction since 1500. We use a near-surface temperature reconstruction on a 0.5°× 0.5°regular grid over [35°-70°N, 25°W -40°E] based on long instrumental series and different proxies (including Greenland ice cores, tree rings and documentary sources) 31,56 . This reconstruction covers the period 1500-2002, although data for 1901-2002 comes from instrumental datasets. Near-surface temperature analyses of the Goddard Institute for Space Studies (GISS, data.giss.nasa.gov/ gistemp/) 32 were herein used at 2°× 2°spatial resolution to update the temperature reconstruction over the period 1901-2019. This monthly observational dataset was herein used at 2°× 2°spatial resolution to update the temperature reconstruction over the period 1901-2019. To do so, reconstructions and instrumental observations were linearly interpolated into a common 2.5°× 2.5°resolution grid over land, as that employed in Barriopedro et al. 57 . Afterwards, seasonal mean temperature anomalies were computed with respect to their respective 1981-2010 climatologies. Finally, the European mean temperature anomaly of each summer in the 1500-2019 period was computed as the area-weighted mean of all land 2.5°grid cells.
To assess whether the decadal frequency of extreme European summers is significantly higher than that expected by random chance we performed a 1000-trial bootstrap, each containing a randomly resampled series of the European summer temperature anomalies over 1500-2019. For each trial, the maximum running decadal frequency was retained, with the 95th percentile of the resulting distribution identifying the value whose one-tailed likelihood of occurring by chance is less than 5%.
C3S soil moisture dataset. The C3S dataset from Copernicus provides estimates of volumetric soil moisture (in m 3 m −3 ) in a layer of 2 to 5 cm depth, retrieved from a large set of satellite sensors. Data is presented on a 0.25°× 0.25°regular grid with some gaps in space and time. Climate Data Records (CDR) and interim-CDR (ICDR) products are generated using the same software and algorithms. CDR is intended to have sufficient length, consistency, and continuity to characterize climate variability and change. ICDR provides a short-delay access to current data where consistency with the CDR baseline is expected but has not been extensively checked. The dataset contains the following products: "active", "passive" and "combined". The "active" and "passive" products are created by using scatterometer and radiometer soil moisture products, respectively. The "combined" product results from a blend based on the two previous products. Here we used the "combined" dataset, which is available for Europe since 1978. Climatological means for each calendar day and grid cell were computed in order to derive local and regional anomalies during the 2019 HW events.
HW algorithm. To perform a spatio-temporal tracking of the 2019 summer HWs, we have adopted a semi-Lagrangian perspective. The 850 hPa temperature (T850) from the NCEP/NCAR dataset was used to analyze the spatio-temporal evolution of extreme temperature patterns, instead of considering HWs as isolated local surface extremes, thus enabling the temporal monitoring of the spatial extent of HWs affecting distinct areas during their lifecycle. The algorithm identifies HW events, defined as areas larger than 500,000 km 2 with daily mean T850 above the local daily 95th percentile (with respect to 1981-2010) that persist for at least four consecutive days and fulfil some predefined conditions on spatial overlap during those days. Additional information on this methodology can be found in Sánchez-Benítez et al. 58 .
Saharan intrusions. A catalogue of air masses with subtropical desertic characteristics was obtained relying on simple thermodynamic air properties, considering the following conditions: (i) 1000-500 hPa geopotential height thickness higher than 5800 m 59 ; (ii) 925-700 hPa potential temperature (θ) above 40°C.
Grid cells satisfying both criteria correspond to low density, warm, stable and very dry air masses, with the potential to be additionally warmed by subsidence 60 . Using the NCEP/NCAR reanalysis dataset, we have classified the mean climatological (1948-2019) location and extension of Saharan air masses during summer, identifying temporary intrusions towards higher latitudes for each grid cell on a daily basis. Further details on the methodology can be found in Sousa et al. 12 .
Temperature tendency and related processes. The contributions of horizontal advection and vertical descent to temperature tendency were determined as where (1) is the temperature advection by the horizontal wind, and (2) the temperature tendency by vertical motion. Equations (1) and (2) are computed from daily mean fields in constant pressure coordinates, according to the pressure levels available in the NCEP/NCAR dataset, with (λ, ϕ, t) representing latitude, longitude and time, respectively, and v being the horizontal wind, T the temperature, ω the vertical velocity and θ the potential temperature. The daily mean temperature rate due to other diabatic processes (e.g. radiative and heat fluxes) is estimated as a residual from the previous two terms based on the temperature tendency equation: where the first term on the right-hand side of (3) is the daily mean temperature tendency (in°C day −1 ). It must be kept in mind that different factors such as sub-grid turbulent mixing, analysis increments and other numerical errors may contribute to the residual term. This bulk analysis is performed for the 1000-850 hPa layer. The relative contribution of each term to the temperature tendency is used to identify the dominant mechanism for each day and grid cell. For further details on the methodology the reader is referred to Sousa et al. 26 .
Analogue method. We use the analogue method, which infers the probability distribution of a target field from the atmospheric circulation during a considered time interval 28 . Herein, two analogue exercises were designed, one for each HW event, but with different target fields, as explained below. In both cases, flow analogue days are defined from their root-mean-square errors (RMSE) with respect to the actual Z500 anomaly field at the time of the HW event over a given domain ([35°-65°N, 10°W-25°E] for the June HW, and [40°-70°N, 10°W-25°E] for the July HW, following the regions with the largest Z500 anomalies; Fig. 2b For each day, the random selection of flow analogues was repeated 5000 times to derive flow-conditioned distributions. To test whether the dynamics played a significant role in the reconstructed anomalies of the target field, unconditional distributions were also retrieved by repeating the whole process with a random selection of days (instead of restricting the search to N days with similar flow configurations). Different choices in the spatial domain or the number of circulation analogues (e.g. N values ranging between 10 and 50) were tested, yielding similar results. For both analogue exercises, Z500 and volumetric soil moisture content at 0-10 cm were obtained from daily means of the NCEP/ NCAR reanalysis (1950-2019). EOBS (1950-2019) was employed for daily maximum temperature at 2 m in the second analogue exercise, although the conclusions remain unchanged if NCEP/ NCAR is used instead. We also repeated the analogue exercises with equivalent fields from different periods and reanalyses: ERA5 (1979-2019) and ERA20C (1900-2010; herein taking the 2019 fields from ERA5). Note that in both ERA reanalyses the uppermost soil layer spans 0-7 cm, and that ERA20C only assimilates surface and mean sea level pressure and surface marine winds. In all cases, anomalies are defined with respect to the 1981-2010 period.
In the first analogue exercise, we reconstructed the expected mean volumetric soil moisture fraction for the 15-day period ( [1,15] day interval) after each day of the June HW (24 June-1 July 2019) by using daily flow analogues from the present and past subperiods separately, defined as 1984-2018 (1999-2018 and 1951-2010) and 1950-1983 (1979-1998 and 1900-1950) in NCEP/NCAR (ERA5 and ERA20C), respectively. This 15-day period is similar to the temporal interval between the end of the June HW and the beginning of the July HW (Fig. 4a), but we obtain similar results for other choices (e.g. [5,20] or [1,30] day intervals). In addition to spatial fields, flow-conditioned and random distributions of the mean soil moisture content over WEU ([43°-53°N, 0°-10°E]) were computed for each subperiod. As the atmospheric circulation is constrained, the difference between the reconstructions of the past and present should largely be ascribed to overall climatological differences between the two subperiods, enabling the estimation of the effect of recent changes in the soil moisture distributions, howsoever caused.
A second flow analogue exercise was performed to address whether the previously accumulated soil moisture deficits over WEU could have contributed to intensifying the temperature anomalies over that region at the time of the July HW. In this case, we reconstructed the maximum 2 m temperature anomalies expected from the circulation during the July HW, distinguishing between analogue days preceded by dry and wet conditions. Wet and dry conditions are defined as summer days of the full period with 15-day mean regional anomalies for the previous [−15,−1] day interval staying above the 66.6 th percentile and below the 33.3rd percentile of the climatological distribution, respectively. That way, soil moisture departures of a given analogue day represent previously accumulated values and are not the direct response to the actual atmospheric circulation conditions. We tested the robustness of the results with respect to the percentile employed for the classification of dry and wet years (e.g. the first and last deciles or quartiles), reporting larger differences for the most extreme definitions. To avoid the effects of long-term trends that may further complicate the causality of the relationships between soil moisture and temperature, these fields were detrended by removing the local trends. For Z500, we removed the regional mean linear trend over the considered domain in order to keep the spatial gradients when searching for flow analogues. Flow-conditioned and random distributions of regional mean temperature anomalies over WEU were also derived for wet and dry conditions. The choice of reanalysis, periods and other methodological aspects can affect the results quantitatively, as well as the spatial details of the reconstructed patterns. However, for the datasets employed and sensitivity tests described above, we did not report substantial differences that affect the main conclusions of the text, therefore adding confidence on the results.