Western Europe’s extreme July 2019 heatwave in a warmer world

Summertime heatwaves are extreme events with a large societal impact. Intensity, duration and spatial extent, all heatwave properties are projected to increase in a warming world, implying that summers that qualified as extreme in the past will become increasingly normal. In this paper we quantify how the changes play out for the July 2019 European heatwave that shattered temperature records throughout Western Europe. We combine a storyline approach with ensemble Pseudo Global Warming (PGW) and high-resolution dynamical downscaling. The downscaling is done with a regional climate model (RACMO2, 12 km resolution) and a convection-permitting model (HCLIM-AROME, 2.5 km resolution). Under PGW the maximum temperature during the heatwave rises 1.5 to 2.5 times faster than the global mean, implying that even at moderate warming levels the heatwave impact changes are tangible. Moreover, there is no sign that the increase in the maximum temperature levels off at higher warming levels, implying that at +4K above present-day temperatures could reach 50 ∘C. During heatwaves cities become islands of heat where daily maxima and night-time minima are up to 5 ∘C higher than in rural areas as we show in ultra-high resolution HCLIM-AROME simulations at 150 m resolution.


Introduction
Historic weather and climate extremes such as summertime heatwaves, winter cold spells, flooding and extreme precipitation events, provide important benchmarking events to society.If impacts were large they can trigger government policy action and their storylines tend to stick in the people's minds for a long time.Impact-modelling applications often require historic observed time series and extremes as input.By transforming the observed time series to a future climate the applications can be stress tested.This so-called 'futurising' can be particularly revealing as these time series provide a means to assess the impact of observed events under future climate conditions.National climate centers often derive these using statistical methods in combination with regional dynamical downscaling (e.g.Bakker 2015, Klein Tank et al 2015, Feigenwinter et al 2018, Bessembinder et al 2023).In Klein Tank et al (2015) for example, a recordsetting rainfall event in the Netherlands was discussed and it was shown how the event could become both more intense and spatially more extensive in a warmer climate (Attema et al 2014, Lenderink andAttema 2015).
The present paper discusses another landmarking event, the European heatwave of July 2019 during which temperature records were broken throughout Western Europe (figure 1).Temperature maxima (Tx) reached values of more than 15 • C above the 1991-2020 climatology (Ma et al 2020).In the Netherlands the 40 • C level was reached for the first time ever, shattering the almost eighty year old national temperature record by more than 2 • .In terms of the wet-bulb globe temperature (WBGT, Liljegren et al 2008) values near the upper limit of level-4 ('High Threat') were reached in the 5-scale category used by the US National Weather Service (see e.g.NWS).Due to the low relative humidity, WBGT did not reach the most severe 'Extreme Threat' level-5 in which a strong recommendation is made to delay all outdoor (work) activities.The July heatwave came after one in June, and barely a month later a third heatwave passed.After the very dry summer 2018 (Aalbers et al 2023), summer 2019 wrote history by its extremely high temperatures.Temperatures that would have been next to impossible without global warming (Vautard et al 2020).
Because of the large amounts by which existing records were broken, the July heatwave triggered research on the maximum possible temperature in the current climate (Fischer et al 2023, Noyelle et al 2023, Zhang and Boos 2023).Where Fischer et al (2023) used ensemble boosting to find whether the occurred situation could have been more extreme, Zhang and Boos (2023) took a more fundamental approach by testing the hypothesis that the maximum attainable temperature is eventually limited by convective instability between surface and mid-troposphere.Noyelle et al (2023) combined this within an analogue framework to unravel the role of thermodynamic and other processes.Both studies concluded that the 50 • C limit cannot be excluded even in the current climate.
In a steady climate the probability of breaking records reduces asymptotically with time.In a warming world however this probability remains non-zero or may even increase (Sterl et al 2008, Rahmstorf and Coumou 2011, Fischer et al 2021).Shattering probabilities can also increase if the relevant underlying distribution changes in higher-order moments.This is seen in the observations as the local summer temperature distribution not only shifts but also changes shape (e.g.Holmes et al 2016, van der Wiel andBintanja 2021).Due to various mechanisms (e.g. the lagged warming of the ocean, changes in the large-scale circulation and vertical stability, enhanced drying of soil moisture and reduced aerosol dimming due to air quality improvements) the hot summer days over much of the European land area warm more rapidly than the global mean and than the local average summer temperature (Haarsma et al 2009, van Oldenborgh et al 2009, Coumou et al 2018).Climate models generally underestimate this factor but do confirm that the extremes increase faster than the mean (Min et al 2013, Vautard et al 2023, Schumacher et al 2024).
Given the previous considerations on the changes in extreme temperatures, it is of interest to find out how the 2019 summer season would change if it occurred in a warmer world.This is the topic of this paper.The following questions are addressed: How different would the July heatwave look and how much higher would maximum temperatures reach?How would the summer season as a whole transform and what are the consequences in terms of cumulative impact measures?How many summers in a warmer climate would start looking similar to 2019?Would 2019 still qualify as exceptional in a warmer world, and at what warming level would it become the norm?
Regional climate models (RCMs) offer the possibility to examine these questions, because we can re-simulate historical extreme weather events, and transform them to future climate conditions using the Pseudo Global Warming (PGW) approach (Schär et al 1996, Brogli et al 2023).In PGW the large-scale day-to-day weather variability is taken from reanalyses data such as ERA5 (Hersbach et al 2020).At the boundaries of the RCM (and to the initial condition) a perturbation-referred to as the delta field-is added.This delta field is typically derived from the future response of a climate model or climate model ensemble such as CMIP6 and aims to include the mean effect of a warmer climate.During the simulation the regional model adapts to the modified conditions, thereby capturing a considerable part of the thermodynamic warming response (Schär et al 1996, Lenderink et al 2019).Because the large-scale circulation conditions are similar in the control and the future simulation, PGW-derived responses generally have a high signal-to-noise ratio (Brogli et al 2019a, 2019b, 2023 2022), the PGW realisations will in terms of circulation be structurally much more similar to their original than the best analogue will likely ever get.Compared to purely statistical methods the PGW dynamic modelling approach has the advantage that all output variables are dynamically consistent.Because we consider extremes the use of an ensemble approach is critical and provides essential information on the uncertainty of the PGW approach.A further motivation for this is that although a heatwave is steered by the large-scale circulation conditions, the actual temperature field and also the future response remain sensitive to small differences in the large-scale forcing (Sousa et al 2018, 2020, Suarez-Gutierrez et al 2020).
In this paper we use the hydrostatic RCM RACMO2 (van Meijgaard et al 2012) for dynamical downscaling the European climate at a horizontal resolution of 12 km, as well as with two recent versions of the non-hydrostatic convection-permitting regional climate model (CPM) HCLIM-AROME (Belusic et al 2020) at 2.5 horizontal resolution.HCLIM-AROME (or simply HCLIM from here) is the climate version of the Harmonie-AROME model used for operational NWP at a number of European weather centers (Bengtsson et al 2017).The two versions of HCLIM used here are HCLIM38 and HCLIM43, the number referring to the underlying model cycle of Harmonie-AROME.As a feasibility study to show to potential of providing very high resolution information to be used for climate change impact assessment, we performed some additional runs at very high resolution (up to 150 m).At the highest resolution one can start to investigate how local warming depends on urbanization, as well as other land surface characteristics.Much of the literature on CPMs initially focused on (extreme) precipitation (Coppola et al 2020, Ban et al 2021, Pichelli et al 2021), but gradually more parameters are being studied, including temperature extremes and heatwaves (Lucas-Picher et al 2021, Sangelantoni et al 2023).By using CPM-modelling at increasingly high spatial resolution city-effects are starting to get resolved.
The outline of the paper is as follows.Section 2 discusses data and methods.Sections 3 and 4 describe the main results and section 5 provides a summary and further discussion.

Observations, target area, heatwave definition
For observations of daily mean (Tg) and daily maximum (Tx) temperature and precipitation, we use the E-OBS data (Cornes et al 2018) at 0.1 degree resolution, as well as observations from meteorological weather stations in the Netherlands.Wind at 850 hPa is taken from ERA5 (Hersbach et al 2020).The key area of interest is indicated in figure 1 and encompasses the south-eastern parts of the Netherlands and a small part of labelled SENL).This region experienced the highest temperatures in the Netherlands during the July heatwave.Definitions of heatwaves differ between countries.We use the Dutch definition of a regional heatwave: a series of at least five days with maximum temperatures above 25 • C of which at least three reach 30 • C. Following this definition the heatwave lasted 6 days (from 22-27 July 2019) when computed from Tx averaged over SENL.

Dynamical downscaling
Dynamical downscaling at a horizontal resolution of 12 km is performed with the RCM RACMO2 (van Meijgaard et al 2008(van Meijgaard et al , 2012) ) (RACMO from here).Boundary and initial-conditions are taken mostly from ERA5 reanalysis (Hersbach et al 2020).As explained below and in appendix A, the control simulations are then transformed to a warmer climate using the PGW method (section 2.3).A number of simulations has been downscaled further using two versions of the non-hydrostatic CPM HCLIM-AROME (Belusic et al 2020).HCLIM-AROME is the climate version of the Harmonie-AROME model used for NWP at a number of European weather centers (Bengtsson et al 2017).Two versions were used, HCLIM38 and HCLIM43, the number indicating the model cycle of Harmonie-AROME.Collectively we will refer to them as HCLIM.The default horizontal resolution of HCLIM is 2.5 km.
As a feasibility study to show the potential of providing very high resolution information to be used for climate change impact assessment, we performed some additional runs with HCLIM43 at very high resolution up to 150 m.In terms of model settings, an important difference is that at the finest resolution the shallow-convection parameterisation (de Rooy et al 2022) was switched off.Presently it is only feasible to run the 150 m configuration on small domains and for short periods.HCLIM in the default setup (2.5 km) is run nested in RACMO2.Soil moisture initialization is then taken from the particular RACMO2 simulation in which it is nested.The HCLIM43 simulations at higher resolution are subsequently nested in the 2.5 km HCLIM43 run (More details in appendix B).

Construction of the PGW ensemble
To transform the historic weather to future climate conditions we use the PGW approach (Schär et al 1996).In this dynamic modelling approach ERA5 provides boundary and initial conditions for the reference control (CTL) simulations of the RCM.In a second simulation (the PGW-simulation) the boundaries and the initial condition of the RCM are perturbed by adding a delta field.This delta field is derived from the future climate response of a single GCM or a multi-GCM ensemble like CMIP6.We explored different delta fields and a range of global warming levels (dT g ) between +1.5K and +4.0K (above present-day, since they are added on the CTL simulation).They are listed in table 1 and referred to as the RACMO PGW-ensemble.The first two groups are based on single-model ensemble-means (Aalbers et al 2023).For the CMIP6-based deltas we included the same 33 models (table A1) that were used in van der Wiel et al (2024) as a basis for constructing a new set of storyline climate scenarios for the Netherlands.The dry11 and wet11 subgroups are subsets with respectively an above-average and below-average future annual precipitation response (appendix A).Finally, two simulations were performed with modified soil parameters (see section 2.4).
The PGW ensemble is based on delta fields derived from GCM responses in the high-end scenarios of CMIP5 (rcp8.5)and CMIP6 (ssp5-8.5).However, by using global warming levels rather than fixed time horizons the simulations are not extreme by default.Moreover, for a given dT g the delta-fields of different SSP scenarios are often quite similar, i.e. there is approximate path independence.We therefore expect our results to be relatively robust to the choice of scenario.Note that in lower emission scenarios it may be difficult or even impossible to reach high values of dT g .

Structure of the delta field
The delta field is constructed for surface pressure, surface temperature and for 3D temperature, humidity and horizontal wind at 19 pressure levels, which are interpolated to the model grid.The monthly-based delta field is interpolated linearly in time to provide a smoothly varying perturbation.Examples of the summer (JJA) are given in figure 2. Prominent common features are the enhanced warming over land and towards the Mediterranean, and a lagged warming of the North Atlantic, with an associated high-pressure response with a center west of Ireland (Haarsma et al 2009(Haarsma et al , 2015)).Most summer delta fields involve a weak easterly circulation component and a stabilisation of the atmosphere, but details differ.The circulation changes are weaker in CMIP6-all33 (top-right panel) than in EC-Earth (top-left) among others because the latter is based on averaging over multiple GCMs.The 'dry-11' subset has a more eastward extension of the high and a stronger circulation response, while the 'wet-11' subset is clearly cooler with a weaker circulation response.An increase in the (vertical) stability of the atmosphere is also seen in many delta fields (not shown).This increased stability comes with accompanying reductions in cloud cover and near-surface relative humidity, and an increase in shortwave radiation.However, note that except from initialisation, the delta fields are only added at the lateral boundaries of the RCM.Note that the actual difference between a future PGW and a control simulation will generally differ from figure 2.

Additional RACMO simulations with a modified surface parameters
RACMO has a systematic low bias in maximum temperature Tx.To investigate whether the bias has a large effect on the Tx future response, we performed additional sensitivity runs in which the surface scheme in RACMO was modified by reducing the root depth by 20% and increasing rapid runoff by 25%.In the summer these modifications have as a result that soil drying limits evaporation more strongly.This leads to a higher Tx and hence reduces the bias.Note that the modifications also influence other parameters like precipitation.In the winter these adjustments have little impact.
The sensitivity experiment serves a second purpose, to examine dependence of the July 2019 heatwave to the soil moisture state.Soil moisture conditions over SENL at the start of the heatwave were low, but not exceptionally low.Due to the more rapid runoff in the modified run, the soil moisture state at the start of the July 2019 heatwave is much drier.The difference between the modified CTL simulation and the original CTL (with default surface settings) is an estimate of the effect of the soil moisture on the intensity of the heatwave.

Bias adjustment
The aforementioned negative bias in Tx necessitates bias adjustment.When discussing future RACMO PGW time-series we apply a simple time-dependent correction and subtract from the PGW simulations the difference between the ensemble-mean CTL and the observations.For the ensemble-mean this is equivalent to adding the future response (PGW-CTL) to the observations.Quantile bias correction (Cannon et al 2015) gave similar results (not shown).

A PGW-based method to transform the observations to a warmer climate
To obtain a future transformed time series of Tx for arbitrary global warming level dT g .First we scale the (time-dependent) PGW-CTL response of each simulation by its dT g value.The result is referred to as the Temperature Response Factor (TRF).The TRF is then added to the observations after multiplication by dT g (essentially a form of pattern scaling) to produce a futurized series x fut (t) = x obs (t) + dT g • TRF(t).Ensemble-spread in the TRF can be interpreted as an estimate of the (instantaneous, flow-dependent) uncertainty of the future projection.At the daily time scale the TRF is relatively noisy due to for example subtle timing differences etc, but at the 3-day timescale shown here this is less of an issue.The above method is simple and relies on the near-synchronous evolution of the RCM and observations.Various other approaches are possible (e.g.Cannon et al 2015, Whan et al 2021, Dieng et al 2022).In quantile-delta mapping (QDM) for example one first determines a number of quantiles of the observed and model distributions.In the second step the future response of each quantile is added to the same quantile in the observations.While QDM has as an advantage that it does not require synchronized time evolution, it cannot account for flow-dependent response uncertainties like the TRF-based method.

Seasonal impact metrics
Three impact metrics are considered.The first is the cumulative heatwave intensity index (CHWI), which is the cumulative temperature excess above 25 • C during regional heatwaves (heatwave definition in section 2.1).For persistence we simply use the total number of days spent in heatwaves (CHWD).Finally we introduce the number of 'cool-off ' days (Tx below 25 • C).To estimate uncertainty in these metrics we use a simple bootstrap method.A bootstrap ensemble was created by randomly drawing at each time instant a TRF value in between the 90% range of available TRF values (at that time) assuming a uniform distribution.For each of these synthetic runs the metrics were computed and a confidence interval was defined from the 90% range of these values.

Warming of the 2019 summer in the PGW-ensemble
The following two sections describe our results.We center discussion around 3-day maximum temperature (Tx3), but will also briefly discuss 3-day mean temperature (Tg3).Results in this section are based primarily on the RCM simulations.The CPM results follow in section 4.

Current climate
Figure 3 shows Tx3 in the SENL region from June to mid September in the observations (black line) and the CTL simulations of the PGW-ensemble.The 3-day precipitation sum is added at the bottom of the figure.A couple of robust features emerge.By design the CTL simulations follow the observations closely in terms of the temperature variability.While ensemble spread varies in time, the ensemble is rather consistent during the peak of the July heatwave.Precipitation is clearly correlated with the ensemble spread in Tx3; whenever Tx3 shows large spread there is usually also considerable spread in the precipitation.RACMO clearly has a cold bias in Tx.This was the primary motivation for the additional RACMO run with modified surface parameters that limit evapotranspiration.This modification (dashed line, label 'Rootdc') has a smaller bias in Tx3 especially during the heatwave(s).Therefore, (excess) evaporation augmented by or following from (too) high soil-moisture availability is one of the contributing factors to the cold bias in Tx.

Soil moisture
The sensitivity experiment also shows that had the soil moisture been lower (and thus evaporation been even more constrained) Tx3 could have been higher still (see also discussion in Noyelle et al 2023).For example around the same date in 2018 the soil moisture was clearly drier.Figure 4 illustrates this.Here Tx3 is shown as a function of the soil wetness index of the topmost meter (swi1m).The colored lines indicate several swi1m-conditional Tx3 percentiles, obtained by binning the data in the x-direction, independently for the original and modified simulation.The percentile lines show that the sensitivity of Tx3 to the value of swi1m in the original and modified simulation are quite similar.Soil-moisture sensitivity of Tx3 increases at the low end of swi1m.The difference of ∼4 • C between Tx3 at the warmest day in July 2019 (colored dots) is higher than the prediction of the highest displayed percentile at that value of swi1m.

Future climate
Figure 5 shows Tx3 in the PGW-ensemble for dT g = 3 K, for which we have most simulations and the runs with amodified surface scheme (dashed line).Figure S2 shows the other warming levels.We subtracted the 3-day difference between the ensemble-mean of the reference simulation and the observations (see Methods).Not surprising the future simulations are nearly always warmer than the observations (An exception occurs around 10 June when the ensemble shows large discrepancy in precipitation).In the warmer climate Tx3 peaks at values above 44 • C during the July heatwave.Locally and at individual days Tx reaches higher values.The TRF (future response divided by the applied warming level, see Methods) reaches values between 1.5 and 2.5 during the heatwaves.
Ensemble variations arise from initial-condition uncertainty (like in the control simulation generally quite small), but more importantly from differences in the delta fields (Aalbers et al 2023).Especially towards the end of the July heatwave the ensemble spread is very large.The explanation is straightforward.The end of the July heatwave was caused by cooler air approaching from the west.Members that use a PGW delta field  with a relatively strong easterly wind component (figure 2) show a delayed end date, members with a weaker circulation component do not.The time-varying spread of the ensemble illustrates that using ensembles is essential when applying PGW to case studies.Note further that the simulation with the modified surface scheme (dashed line) generally has an above-average response in early to mid summer, but is not always an outlier.Results are not very sensitive to the size or location of the target region (figure S3 illustrates this for a region in France).

Temperature response factor (TRF) during the July heatwave
Figure 6 shows the average Tx response during 22-27 July 2019, the period of the heatwave in the observations (for definition, see Methods).Lines of constant TRF are indicated.Note that not all models shown in the legend have simulations at all warming levels and also that in order to increase visibility each point has been shifted by a small random amount around the warming level (horizontal axis).As earlier mentioned, initial-condition uncertainty (i.e. the vertical line segments if available) is small (⩽1 degree), compared to the impact of the delta field.At all warming levels the smallest response is seen in the simulations labelled CMIP6-wet11.The delta field of these simulations was based on a selection of 11 CMIP6 models that show a relatively weak drying signal (precipitation decrease) over The Netherlands and  Oblique lines indicate TRF-levels between 1 and 2.5.To increase visibility each point has been shifted by a small random amount on the horizontal axis.The lowest three simulations in the legend form an initial-condition ensemble (11 members per warming level).Their spread (±1.64σ) is indicated by the vertical line segments.Note that not all models shown in the legend have simulations at all warming levels.surrounding countries.In contrast, the TRF of the CMIP6-dry11 selection (based on the 11 CMIP6 models that showed the strongest drying) is among the largest.
The HCLIM simulations at default resolution (2.5 km) are also included.They have a TRF in the same O(1.5-2.5)range as RACMO.As already discussed, the RACMO run with modified surface scheme (labelled RACMO-f-CMIP6-dry11-Rootd) has a Tx response comparable to the equivalent run with the default surface scheme (RACMO-f-CMIP6-dry11).

Mean temperature
Aging societies are particularly vulnerable to high night-time (van Oldenborgh et al 2019, He al 2022).Therefore the analysis has carried for Tg3 figure S4).Observations and the CTL simulations of Tg3 are in better agreement than for Tx3.Tg3 reaches values around 29 • C in the current climate, but this increases to nearly 35 • C under dT g = 3.0 K, implying a TRF of around 2. The RACMO simulation with modified surface parameters has a positive Tg3 bias, but like Tx3 its future response is similar.

Seasonal impact metrics
Figure 7 summarises the three impact metrics aggregated over the summer of 2019 a function of global warming level.Note that here we first applied the TRF-method to derive future time series for arbitrary global warming level (figure S5 confirms that the series obtained in this way are indeed similar to figure 5).The light grey shading shows the envelope of the most optimistic and most pessimistic scenario (obtained by using a low or high value of the TRF at all days).Dark grey shading indicates the 90% range of a simple bootstrap ensemble (see Methods).The width of the envelope (i.e. the uncertainty) increases with warming level because it is computed from the spread of the TRF multiplied by the warming level.
There is a gradual increase in the intensity and duration and a similar decrease in number of cool days as the climate warms.At the highest warming levels the summer of 2019 becomes a near-continuous heatwave with very little time to recover in between the heatwaves.Already at a +2K warming level the cumulative intensity index (left) increases by at least 40% (P05 lower bound i.e. the most optimistic scenario), but on average more than a factor 2. Despite the hot summer there were still many more cool days (Tx ⩽ 25 • C) than days spent as part of a heatwave.This cool-off time shortens progressively and at +2K global warming level the number of heatwave and cool-off days become approximately equal.Similar tendencies are found (but with different numbers) for other summers (we return to this point in the discussion).

Convection-permitting modelling
Up to now we mostly discussed statistics for the target region.These statistics were mostly derived from the RCM simulations.Now we look at the spatial pattern, for which we use the CPM data.As explained in appendix B we use a HCLIM43 simulation that was initiated on 1 July 2019 using as boundary forcing the PGW+2K RACMO output from the CMIP6-based dry11 subset.).The main argument for using this 'late-start' simulation is that performed well in terms of Tx (results below).

Spatial pattern and heatwave area
Figure 8 shows the time evolution and the spatial pattern.Unlike RACMO (which was generally too cold in maximum temperature) HCLIM43 attains maximum temperatures close to the observations over the SENL region (Averaged over the period 14-31 July the bias in Tx over SENL is just 0.04 • C; for completeness this value has been subtracted from both CTL and PGW simulations using the same value for the entire domain.This seems appropriate for a region around the Netherlands).
With Tx3 above 42.5 • C in the PGW+2K simulation, the associated TRF ⩾ 2 over SENL (panel (b)).The bottom row show the spatial pattern of the 3 warmest days (24-26 July).For E-OBS the highest values are found along the Rhine valley just across the German border.The HCLIM43 result confirms these higher values.Under PGW+2K (panel (f)) almost the entire region becomes warmer than the warmest locations in the current climate.A huge area experiences a Tx3 (well) above 40 • C (this area increases almost exponentially with dT g , see also the expanded figure S6).Averaged over the 6 days of the heatwave the TRF (panel (c)) was above 1.5 in many locations in Western Europe, reaching values up to 2.5 in southern parts of the Netherlands.Associated with the high pressure response west of The Netherlands in the PGW forcing pattern, cooler air reaches the northern parts of the Netherlands and Germany, causing TRF values around 1 in those areas.

Urbanization
At the scale of Western Europe, the spatial patterns of E-OBS and HCLIM43 look quite similar (figures 8(d)-(f)).However, at a smaller scale HCLIM43 displays much more spatial inhomogeneity.The darkest spots (i.e.highest temperatures) typically occur in regions with a high degree of urbanization.At 2.5 km spatial resolution regional climate modelling with HCLIM is clearly starting to resolve city effects.Contrast this to the E-OBS gridded data set.E-OBS is generated by spatially interpolating data from the meteorological stations (circles in the maps) to a 0.1 • gridded product.The meteorological stations are typically positioned outside the main cities.As a result, the cities do not really show up in E-OBS (except along the German sector of the Rhine).On the other hand, it is well known that cities experience higher maxima than the rural areas around it.HCLIM starts resolving city-effects at 2.5 km resolution, but at even higher spatial resolution, the difference between rural and urban areas becomes much more visible, as discussed next.

Convection-permitting simulations at 150 m resolution
Finally we look at the 150 m resolution HCLIM43 simulations.On a small domain it is even more essential to provide the right boundary conditions.For this reason the 150 m runs are nested in the 'late-start' HCLIM43 simulations (details in appendix B).Like for the other HCLIM maps (figure 8) we subtract a single bias value, namely the mean difference wrt E-OBS averaged over (the part of) SENL covered by HCLIM43 over the 150 m simulation period 21-31 July 2019.For Tx the bias is 0.6 • C. It needs to be emphasised that the 150 m resolution simulations are experimental.Knowledge HCLIM in this configuration is limited.The output is used to illustrate typical differences between urban and rural areas.
Figure 9 shows Tx averaged over the three hottest days.Although qualitatively similar to the 2.5 km result, the imprint of the underlying land type is much stronger.A black contour is added where the urban fraction is 50%.In urban and industrialised areas (e.g. in the western part of the domain) the maximum temperatures are typically 1-5 • warmer than in nearby more rural areas.Areas with a sand type underground are also clearly warmer than the low-lying polders in the west.HCLIM locally reaches higher temperatures than E-OBS and even in CTL Tx3 reaches 40 • C in many places.The right panel shows the  PGW +2K result.Like figure 8(a) large area experiences Tx3 above 45 • C, confirming a response much larger than both dT g and the mean delta field.
In the evening and during the early hours of the night cities remain much warmer than the rural areas especially during heatwaves.At night the cities finally cool down, but remain warmer than the neighbouring areas.Supplementary figure S7 shows Tn3, the 3 day average minimum temperature during 24-26 July.Large regions show Tn3 well over 20 • C. In the PGW simulation a Tn3 of 25 • C is no longer exceptional (figure S8).Such tropical nights lead to increased mortality rates, as they interrupt the normal sleep physiology (He et al 2022).
One of the primary motivations for doing the high-resolution simulations is to find out whether cities experience a similar (or weaker or stronger) response than more rural areas.As the response patterns display less spatial heterogeneity than the CTL and PGW patterns themselves (figure S9), this implies that cities not necessarily experience a stronger response than more rural areas.To make this a bit more explicit and quantitative, figure 10 shows histograms of Tx3 during 24-26 July for the different subsets of points.For land points the Tx3 is on average around 2.5 • C warmer in urbanized areas than in rural areas.Within the rural subset (middle panel) a double peak in the distribution hints on systematic differences due to the underlying land types.Here we see even more clearly that the scaled response is on average higher in the rural areas, a point that calls for further research.For Tn3 the response is smaller than for Tx3 (figure S8) and appears to be more evenly spread between urban and rural regions.
Finally, note that in the PGW simulations land-use changes are not taken into account.Cities are assumed to preserve their original sizes, undoubtedly a crude approximation.When cities expand, so will the red spots on the maps (figure 9).Current research in high-resolution climate modelling addresses questions like to what extent green roofs for example could help cities to prevent from overheating.

Summary and discussion
In this paper we discussed how the intense heatwave of July 2019 could look if it were to occur in a warmer climate.This is done using regional climate modelling and the PGW approach.We summarise our conclusions and expand on a couple of points.
Intensity, duration and spatial extent.We focused on three defining characteristics of heatwaves: intensity, duration and spatial extent.All three will change (along with increases in frequency).While the quantitative details will vary from case to case, and from year to year, some aspects of the 2019 results are quite robust.The intensity increases more strongly than the imposed level of global warming: For every degree of global warming the modelled maximum temperature typically increases by 1.5-2.5 • .According to the PGW experiments conducted, there is no real flattening of this rate at the higher global warming levels (figure 6).This implies that even at moderate global warming levels the changes are substantial.For example, at dT g = +2K the heatwave duration over the SENL region increases from 6 to 17 days (figure S5).Ensemble spread caused by initial-condition uncertainty is relatively small compared to differences due to the applied delta field.Systematically higher temperature responses were found for a delta field based on the 'dry-11' CMIP6 subset; these simulations warmed more strongly than those based on the 'wet-11' CMIP6 subset (e.g.figure 6).Heatwave duration increases to first order simply as a direct consequence of the higher temperatures (i.e. an earlier start and a later end).Further persistence changes could occur if the large-scale conditions are supportive, as we saw in several members of the PGW ensemble.However, the PGW approach is not well suited to give definite answers to the question of the duration-change, because the large-scale variability is fixed (and equal to the historic conditions).Areal extent-computed as the average area where the 40 • C threshold is exceeded-increases almost exponentially.Finally we discussed differences between cities and more rural areas, by looking at experimental simulations at ultra-high spatial resolution.Already in the current climate the temperatures are substantially higher in the major cities, especially during heatwaves.This difference will remain in a warmer climate, implying that (large) cities will be the first to face extremely high temperatures, both during daytime and during the nights.
Representativeness of the July 2019 heatwave for other hot days.With its record-shattering temperatures the July 2019 summer heatwave was clearly exceptional in the current climate.Since there is a large variation in the heatwave-parameters not all future heatwaves will reach the same intensities as the futurized July 2019.Are there any general conclusions to be drawn from 2019?One aspect is the scaled temperature response (TRF), the response divided by the (global) warming level.During most of the summer of 2019, it was above unity and (generally) higher on the warmest days.To check the representativity of the July 2019 case for other hot days we computed the scaled response using May-September data from all available years and simulations (table 1). Figure S10 shows boxplots of the TRF of Tx3, binned by the observed Tx3 in the target region.The TRF increases with temperature in both RACMO and HCLIM and for the majority of the warmest observed temperatures exceeds 1.5, similar to the response of the 2019 heatwave.Note the mean temperature in the region also increases faster than the global mean.On cold days the TRF is generally below unity, which is not surprising since these days occur typically when the wind is blowing from the west or north, which is the region where the warming lags behind (figure 2).Similar values were obtained for a quantile-quantile analysis of a conventional RCM ensemble (figure S11).Based on these results over all years and simulations we are confident that heatwaves get warmer at a rate 1.5-2.5 times faster than the rise in global temperature.
Would the present-day 2019 summer still be exceptional in a warmer climate?A first-order answer to this question can be obtained by creating futurized time series of all historically available years (we use 1951-2021) and by comparing these to the historic series.Here we simply add the median scaling factor (horizontal lines in figure S10) to the observed daily time series and compute statistics like heatwave intensity, duration, and cumulative heatwave-number.Figure 11 shows the result for a +2K scenario.In terms of the total summed heatwave-number (CHWI) the summer of 2019 as it occurred ranked second after 2018, out of the 71 summers .In the +2K warmer climate, the 2018 summer as it occurred would no longer appear in the Top-10 in terms of CHWI, and 2019 would be lower still.In terms of maximum temperature, the year 2019 as it occurred would still rank in the top-10% in the +2K warmer climate.Similar changes could also be found for a lower warming level of +1.5K but a higher assumed response (P75).Not surprisingly the results get much more extreme at stronger warming levels (figure S12).But even under scenarios with relatively moderate additional global warming, summers like 2018 and 2019 that are in our memory of having been exceptional become increasingly normal.PGW method and alternatives.To conclude we briefly reflect on the PGW method.An important ingredient is the delta field.This delta field describes the average future response of (ensembles of) GCM simulations such as (subsets of) CMIP6.In the interior the regional model adapts to the modified boundary conditions.By doing this, the main thermodynamic response is included, under the assumption that large-scale variability changes are small (Brogli et al 2023).Although the time correlations between the control and PGW are quite high for temperature an ensemble approach is strongly recommended (Aalbers et al 2023).Recently Sanchez-Benitez et al (2022) proposed an alternative approach to futurization, by spectrally nudging a GCM to the observed ERA5 circulation.Also in their future scenarios the July 2019 heatwave warms strongly and they conclude that temperatures around 50 • C cannot be excluded in Western Europe if July 2019 would occur in a globally +4K warmer world.Our results point in the same direction.In fact, recent studies suggest that such high values can perhaps not even be excluded for the region in the current climate (Noyelle et al 2023, Zhang andBoos 2023).More research is required to assess the limit of applicability of PGW approaches, and to compare it to other conventional or new methods.frequency.Because running HCLIM is computationally demanding, only a small subset of the RACMO simulations has been downscaled.In this study we use: (i) HCLIM38 (PGW+1.5K,rcp8.5):These simulations use ERA-Interim (2008-2017) and ERA5 (2018-2020) as boundary and initial conditions.The delta field was based on EC-Earth2.3and a warming level of +1.5K.The horizontal resolution of HCLIM38 is 2.5 km.HCLIM38 is described in Belusic et al (2020).(ii) HCLIM43 (PGW+2K, ssp5-8.5):These simulations use ERA5 as boundary and initial conditions.The delta field was based on the CMIP6 'dry-11' subset.Runs were initiated on 1 January 2019.A second 'late-start' simulation was initiated on 1 July 2019.The horizontal resolution of the HCLIM43 runs is 2.5 km.The main code is shared with Harmonie Cy40 used for NWP (Bengtsson et al 2017), but in HCLIM43 (and Harmonie Cy43) there are differences in the turbulence, convection and cloud scheme (de Rooy et al 2022), as well as in the microphysics (Contreras Osorio et al 2022) and the applied surface physiographical database (notably in the surface/soil scheme).The main difference between HCLIM43 and NWP Cy43 is the applied surface scheme, namely a force-restore scheme in NWP and a diffusion scheme in HCLIM43.(iii) HCLIM43 ultra high-resolution (PGW+2K, ssp5-8.5):These highest resolution simulations (500 m-150 m horizontal resolution) were run on a smaller domain and initiated on 19 July 2023 (i.e.only a few days prior to the heatwave event) and were nested within the 2.5 km resolution 'late-start' HCLIM43 runs.
The reason for using the 'late-start' HCLIM43 runs in the storyline approach is that it had the smallest bias compared to observations.HCLIM is quite sensitive to soil-moisture drying in spring and early summer during conditions that favour evaporation (e.g.abundant sunshine, and anomalously high temperatures).The most likely explanation for the drying is an overestimation of the latent heat fluxes which might be related to details of physiography-related parameters like leaf-area index.If the soil-moisture is below average and a heatwave happens to come by, HCLIM has a tendency to produce too high maximum temperatures.The physical reason is similar to what we saw in the RACMO simulation with the modified soil scheme: reduced soil moisture reduces evaporative cooling and thereby promotes higher temperatures.In the late-start HCLIM43 simulations the soil-moisture state is closer to neutral/climatology, and the bias in the maximum temperature during the heatwave is small.Note incidentally that the HCLIM43 simulation that was initiated on 1 January 2019 had lower soil-moisture at the start of the July heatwave and developed a positive bias in Tx3 (not shown), while displaying a smaller response (figure 6).A more complete analysis of the temperature-bias in terms of the underlying processes is reserved for another study.

Figure 1 .
Figure 1.(left) Maximum temperature (Tx) from E-OBS on 25 July 2019 (shading; the 40 • C isotherm is indicated by a light contour).Streamlines show mean wind at 850 hPa (ERA5).The dashed box is the target region (SENL).(right) Amount of record-exceedance in Tx in 2019 (not all occurred during the July 2019 heatwave).

Figure 3 .
Figure 3. Tx3 in SENL in observations (black line) and in RCM control simulations.Also shown is 3-day precipitation sum (right y-axis).

Figure 4 .
Figure 4. Tx3 when above 20 • C as a function of soil wetness index of the topmost meter in SENL, in the original CTL simulations and those with surface-parameter adjustments (labelled 'ctl rootd') for the period 1979-2020.The lines indicate various conditional percentiles of Tx3, derived from binning in the x-direction.

Figure 5 .
Figure 5.As in figure 3 but for the PGW+3K model simulations.A bias adjustment correction ('BC' in the top-left legend) has been applied to Tx3 (see methods).

Figure 6 .
Figure 6.Response in Tx over the SENL region during 22-27 July 2019, as a function of global warming level.Oblique lines indicate TRF-levels between 1 and 2.5.To increase visibility each point has been shifted by a small random amount on the horizontal axis.The lowest three simulations in the legend form an initial-condition ensemble (11 members per warming level).Their spread (±1.64σ) is indicated by the vertical line segments.Note that not all models shown in the legend have simulations at all warming levels.

Figure 7 .
Figure 7. Cumulative changes over the summer of 2019 (Period covered: 1 June-15 September), as a function of the global warming level.Left: cumulative heatwave intensity; Middle: cumulative number of days spent in heatwaves; Right: number of cool days with daily max temperature below 25 • C. Light shading indicates most optimistic/pessimistic scenarios, dark shading the variations obtained from a bootstrap ensemble (see text).

Figure 9 .
Figure 9. Average Tx during 24-26 July 2019 in the 150 m resolution HCLIM43 simulation.colored circles indicate synops-stations.A thin black contour is used where model urban fraction is 50%.Left: CTL simulation.The black line is the simulation domain; E-OBS is shown outside this box.Right: PGW+2K simulation.

Figure 11 .
Figure 11.Heatwave diagram for CTL and PGW+2K warming scenario assuming a median response.Black/red circles show years with at least one heatwave in CTL/+2K climate; the radius equals the May-September CHWI scaled by the historic mean.Shaded boxes show the 90% percentile (P90) in each direction.Years are labeled if and of the three metrics exceeds its P90.The top-right table shows the top-10 CHWI-ranked years in CTL.The bottom-right table shows the years in PGW with a higher CHWI than the historic maximum (filled circles indicate years absent in the top-right table).

Table 1 .
The RACMO PGW simulations used in this study.
Details: a Simulations forming an initial-condition ensemble (11 members for each warming level and forcing).b RACMO runs with modified soil parameters.