Impact of groundwater representation on heat events in regional climate simulations over Europe

The representation of groundwater is simplified in most regional climate models (RCMs), potentially leading to biases in the simulations. This study introduces a unique dataset from the regional Terrestrial Systems Modelling Platform (TSMP) driven by the Max Planck Institute Earth System Model at Low Resolution (MPI-ESM-LR) boundary conditions in the context of dynamical downscaling of global climate models (GCMs) for climate change studies. TSMP explicitly simulates full 3D soil and groundwater dynamics together with overland flow, including complete water and energy cycles from the bedrock to the top of the atmosphere. By comparing the statistics of heat events, i.e., a series of consecutive days with a near-surface temperature exceeding the 90th percentile of the reference period, from TSMP and those from GCM–RCM simulations with simplified groundwater dynamics from the COordinated Regional Climate Downscaling EXperiment (CORDEX) for the European domain, we aim to improve the understanding of how groundwater representation affects heat events in Europe.

The analysis was carried out using RCM outputs for the summer seasons of 1976–2005 relative to the reference period of 1961–1990. While our results show that TSMP simulates heat events consistently with the CORDEX ensemble, there are some systematic differences that we attribute to the more realistic representation of groundwater in TSMP. Compared to the CORDEX ensemble, TSMP simulates fewer hot days (i.e., days with a near-surface temperature exceeding the 90th percentile of the reference period) and lower interannual variability and decadal change in the number of hot days on average over Europe. TSMP systematically simulates fewer heat waves (i.e., heat events lasting 6 d or more) compared to the CORDEX ensemble; moreover, they are shorter and less intense. The Iberian Peninsula is particularly sensitive with respect to groundwater. Therefore, incorporating an explicit 3D groundwater representation in RCMs may be a key in reducing biases in simulated duration, intensity, and frequency of heat waves in Europe. The results highlight the importance of hydrological processes for the long-term regional climate simulations and provide indications of possible potential implications for climate change projections.

1 Introduction

Over the past few decades, the number of heat waves has increased (e.g., Frich et al., 2002; Christidis et al., 2015; Zhang et al., 2020). The years 2003, 2010, 2018, 2022, and 2023 were among the hottest in Europe, characterized by record-breaking near-surface air temperatures (e.g., Stott et al., 2004; Barriopedro et al., 2011; Dirmeyer et al., 2021; Yule et al., 2023; Jiang et al., 2024; Lemus-Canovas et al., 2024). With projected climate change, the occurrence of heat waves will continue to increase (e.g., Russo et al., 2015; Hari et al., 2020; Molina et al., 2020; Masson-Delmotte et al., 2021), leading to multiple negative socio-economic impacts (e.g., Amengual et al., 2014; Yin et al., 2022).

The physical mechanisms underlying heat waves have been extensively studied (e.g., Horton et al., 2016; Liu et al., 2020; Barriopedro et al., 2023). Heat waves are triggered by strong, persistent, quasi-stationary large-scale high-pressure systems associated with atmospheric blocking events, resulting in subsiding adiabatically warmed air masses and clear skies, allowing for high insolation (Tomczyk and Bednorz, 2016; Kautz et al., 2022). The evolution of heat waves depends primarily on the synoptic weather patterns in combination with ambient soil moisture conditions, further altered by multiple land–atmosphere feedback processes (e.g., Fischer et al., 2007).

European summer heat waves were often preceded by a precipitation deficit in spring (e.g., Stegehuis et al., 2021; Hartick et al., 2021). Due to the long-term soil moisture memory effect, the lack of precipitation in early spring causes negative soil moisture anomalies in early summer and strengthens the land–atmosphere coupling (a measure of the response of the atmosphere to anomalies in the land surface state) with a lower evaporation fraction. In turn, this reduces latent cooling and amplifies summer temperatures. Note that the soil moisture memory is a phenomenon of persistence of wet or dry anomalies over a long period of time, from weeks to months, after the atmospheric conditions that caused them have passed; this allows us to preserve the hydroclimatic conditions of the preceding months (e.g., Song et al., 2019). Thus, the long-term soil moisture memory can contribute to either buffering negative droughts impacts and weakening a heat wave or, conversely, delaying drought recovery and exacerbating the occurrence of a heat wave (Erdenebat and Tomonori, 2018; Martínez-de la Torre and Miguez-Macho, 2019; Hartick et al., 2021; Dirmeyer et al., 2021). In addition to precipitation, soil moisture is strongly influenced by groundwater dynamics via vertical fluxes across the water table (capillary rise) and via horizontal fluxes through gravity-driven lateral transport within the saturated zone. Here, the water table depth dictates the intensity of shallow groundwater–soil moisture and evaporation coupling (Kollet and Maxwell, 2008).

In the context of climate impact assessments, dynamical downscaling of global climate models (GCMs) with regional climate models (RCMs) is widely used to generate regional climate change scenario information (e.g., Mearns et al., 2015; Jacob et al., 2020). RCMs have been shown to provide added value to driving GCMs by better capturing small-scale processes (Giorgi and Gutowski, 2015; Torma et al., 2015; Prein et al., 2016; Rummukainen, 2016; Iles et al., 2020), but model biases (offset during the historical period against observations) and uncertainties in climate projections still remain (Hawkins and Sutton, 2009; Sørland et al., 2018; Evin et al., 2021). In fact, many RCMs tend to overestimate duration, intensity, and frequency of heat waves (e.g., Vautard et al., 2013; Plavcová and Kyselý, 2016; Lhotka et al., 2018; Furusho-Percot et al., 2022).

The role of soil moisture in modelling heat waves is crucial (Seneviratne et al., 2006, 2010; Fischer et al., 2007), but due to the complexity of the processes involved and related high computational cost, the explicit representation of hydrological processes is oversimplified in most RCMs. Commonly applied hydrology schemes are based on 1D parameterizations in the vertical direction with runoff generation at the land surface and a gravity-driven free-drainage approach as the lower boundary condition. In such a parameterization there is no lateral subsurface flow and only the 1D Richards' equation is solved (e.g., Niu et al., 2007). RCMs with a simplified representation of hydrological processes have difficulty reliably reproducing the partitioning of land surface energy flux and, consequently, near-surface air temperatures, leading to warm biases (e.g., Barlage et al., 2021). Hydrological parameters tuning (e.g., Teuling et al., 2009; Bellprat et al., 2016) or developing new parameterizations of groundwater dynamics (e.g., Liang et al., 2003; Schlemmer et al., 2018) have been shown to improve model results. Feedback mechanisms between groundwater, land surface, and atmosphere are also often simplified in RCMs. A physically consistent description of hydrological processes in RCMs can be achieved by an explicit representation of 3D soil and groundwater hydrodynamics together with overland flow, accounting for feedback loops over the terrestrial system and complete water and energy cycles from groundwater across the land surface to the top of the atmosphere (Maxwell et al., 2007), as for instance in the Terrestrial Systems Modelling Platform (TSMP) (e.g., Shrestha et al., 2014; Gasper et al., 2014), a regional climate system model.

Keune et al. (2016) demonstrated a relationship between the representation of groundwater dynamics and near-surface air temperature for the August 2003 European heat wave from the TSMP simulation nested within the ERA-Interim reanalysis (Dee et al., 2011). In their study, the TSMP model was set up over the European domain of the COordinated Regional Climate Downscaling EXperiment (CORDEX) (e.g.,  Gutowski et al., 2016; Jacob et al., 2020) with two different groundwater configurations: simplified 1D free-drainage approach and 3D physics-based variably saturated groundwater dynamics. The clear impact of groundwater dynamics on the land surface water and energy balance is shown: latent heat fluxes are higher and maximum temperatures are lower, especially in areas with shallow water table depth, in the 3D configuration compared to the simplified 1D free-drainage approach. The work of Keune et al. (2016) suggests that 3D groundwater dynamics in TSMP alleviate the evolution of a single heat wave due to weaker land–atmosphere feedbacks compared to a simplified 1D free-drainage approach, at least during the investigated European heat wave of summer 2003. The ability of an explicit representation of groundwater dynamics to moderate air temperatures during a single seasonal heat wave in RCM simulations was also demonstrated in Barlage et al. (2015, 2021) and Mu et al. (2022).

Further studies were carried out to understand whether the aforementioned effects of the groundwater representation persist over longer time periods in RCM evaluation runs and how this manifests itself for heat waves over Europe. Furusho-Percot et al. (2019) show that TSMP simulations driven by the ERA-Interim reanalysis capture climate system dynamics and the succession of warm and cold seasons at the regional scale for PRUDENCE regions (Christensen and Christensen, 2007) consistent with the E-OBS observations (Cornes et al., 2018) for the investigated period of 1996–2018. Moreover, the TSMP multiannual evaluation run exhibits lower deviations of summer heat wave indices from the E-OBS observations compared to the CORDEX RCMs with a simplified representation of groundwater, which tend to simulate too persistent heat waves (Furusho-Percot et al., 2022). This particular behaviour of TSMP is attributed to its improved hydrology, which leads to a better capacity to sustain soil moisture and therefore a more reliable latent heat flux and evaporation. This leads to a decrease in the number of days with anomalously high near-surface temperatures, as well as the intensity and spatial extent of heat waves. An important question still remains: how will these findings be reflected in long-term regional climate simulations in the context of dynamical downscaling of GCMs by RCMs for climate change studies over Europe?

In this paper, we present a unique dataset from the TSMP regional climate system model driven by the Max Planck Institute Earth System Model at Low Resolution (MPI-ESM-LR) (Giorgetta et al., 2013) historical boundary conditions over the European CORDEX domain. We interrogate the statistics of the characteristics of heat events (duration, intensity, frequency) for the summer seasons of 1976–2005 relative to the reference period 1961–1990, by comparing TSMP results and those of the CORDEX RCMs with a simplified representation of groundwater driven by GCM control simulations of phase five of the Coupled Model Intercomparison Project, CMIP5 (Taylor et al., 2012). We strive to better understand the impact of 3D groundwater dynamics on simulated heat events in historical regional climate simulations and potential consequences for climate change projections. While Furusho-Percot et al. (2022) examined the statistics of heat events in the TSMP evaluation run, the long-term TSMP historical climate simulations forced by MPI-ESM-LR GCM have not been previously presented or analysed. Thus, this is the first assessment of the heat event statistics over Europe from a dynamically downscaled GCM with a fully coupled RCM that comprises an explicit representation of groundwater dynamics and two-way non-linear feedbacks from groundwater across the land surface to the top of the atmosphere.

Section 2 introduces the methods, describing the TSMP modelling platform and its set-up, the procedure for detection and analysis of heat events, and the CORDEX ensemble used for comparison with the TSMP results. In Sect. 3, we examine the new TSMP dataset driven by MPI-ESM-LR GCM for consistency with the CORDEX ensemble and present results on the impact of an explicit groundwater representation on simulated heat events in long-term regional historical climate simulations. Section 4 contains the discussion, and Sect. 5 provides the summary and overall conclusions.

2 Methods

2.1 The TSMP modelling platform

TSMP is a scale-consistent, highly modular, fully integrated soil–vegetation–atmosphere regional climate system model (e.g., Shrestha et al., 2014; Gasper et al., 2014). TSMP consists of three component models: the atmospheric COnsortium for Small Scale Modelling (COSMO) model version 5.01, the Community Land Model (CLM) version 3.5, and the hydrological model ParFlow version 3.2. The component models are externally coupled via the Ocean Atmosphere Sea Ice Soil (OASIS) Model Coupling Toolkit (MCT) version 3.0 (Valcke, 2013), which enables interactions between different compartments of the geoecosystem, explicitly reproducing feedbacks in the hydrological cycle from the bedrock into the atmosphere.

COSMO is a non-hydrostatic limited-area atmospheric model (Baldauf et al., 2011). COSMO is based on the primitive thermo-hydrodynamical Euler equations formulated in rotated geographical coordinates and generalized terrain-following height coordinates, which describe compressible flow in a moist atmosphere. COSMO parameterization schemes cover various physical processes, such as radiation, cloud microphysics, and deep convection. The boundary conditions for COSMO are provided by a coarse-grid model, i.e., reanalysis or GCM, whereas the lower boundary conditions (e.g., surface albedo, energy fluxes, surface temperature, surface humidity) are provided by CLM in the current TSMP configuration.

CLM is a biogeophysical model of the land surface (Oleson et al., 2004, 2008). It simulates land–atmosphere exchanges in response to atmospheric forcings. CLM consist of four components that describe biogeophysics, hydrologic cycle, biogeochemistry, and dynamic vegetation. In TSMP, CLM receives short-wave radiation, wind speeds, barometric pressure, precipitation, near-surface temperature, and specific humidity from COSMO. In turn, CLM sends infiltration and evapotranspiration fluxes of each soil layer to the ParFlow hydrological model.

ParFlow is a hydrological model that simulates variably saturated 3D subsurface hydrodynamics using a Richards equation integrated with shallow overland flow based on a kinematic wave approximation (Maxwell and Miller, 2005; Kollet and Maxwell, 2006; Kuffour et al., 2020). ParFlow allows 3D redistribution of subsurface water in a continuum approach. In the TSMP configuration used, ParFlow replaces the hydrological functionality of CLM.

The evaluation run of TSMP was performed by Furusho-Percot et al. (2019) with atmospheric forcings derived from the ERA-Interim reanalysis, which was validated by comparing temperature and precipitation with E-OBS and column water storage with the Gravity Recovery and Climate Experiment (GRACE) satellite data (Landerer et al., 2020). In the recent publication of Ma et al. (2022), the TSMP water table simulation results were used in a machine learning approach and compared to in situ water table observation anomalies over Europe; the results showed good agreement considering that the TSMP model has not been calibrated.

2.2 The TSMP simulation set-up

The TSMP simulation is conducted over the European domain according to the CORDEX simulation protocol (Gutowski et al., 2016) on a rotated latitude–longitude grid with a horizontal resolution of 0.11° (EUR-11) or about 12.5 km for the historical time period from December 1949 to the end of 2005. CLM and ParFlow are initialized with the moisture conditions on 1 December 2011, extracted from the TSMP evaluation run (Furusho-Percot et al., 2019). The COSMO configuration resembles that of the COSMO model in CLimate Mode (CCLM) (Rockel et al., 2008). COSMO extends vertically up to 22 km, divided into 50 levels. CLM has 10 soil layers with a total depth of 3 m, which coincide with the 10 top layers of ParFlow. In addition, ParFlow has five bedrock layers increasing in thickness towards the bottom of the model domain to a total depth of 57 m. The time step for ParFlow and CLM is 900 s, while for COSMO it is 75 s. The coupling time step between TSMP component models is 900 s. The TSMP output constitutes terrestrial essential climate variables with a time step of 3 hours; for details, see Poshyvailo-Strube et al. (2023). The first 10 years of TSMP simulation are discarded due to hydrodynamic spin-up.

Forcing data for COSMO are provided by the Max Planck Institute's MPI-ESM-LR r1i1p1 CMIP5 GCM with a resolution of T63L47 (Giorgetta et al., 2013). For CLM, plant functional types (PFTs) are taken from the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover dataset (Friedl et al., 2002). Leaf area index, stem area index, and monthly bottom and top heights of each PFT are calculated based on the global CLM surface dataset (Oleson et al., 2008). Topography in ParFlow is represented by slopes estimated from the United States Geological Survey GTOPO30 (Daac, 2004). In this study, we improved the representation of subsurface hydrogeology in ParFlow compared to the previous studies (e.g., Furusho-Percot et al., 2019, 2022), where soil parameters were assumed to be vertically homogenous. Here, an aquifer network is added to ensure the relationship between surface and subsurface water flow (Naz et al., 2023). Static land surface input data, which includes soil colour, percentage of clay and sand, dominant land use type, dominant soil types in the top layers, dominant soil types in the bottom layers, and subsurface aquifer and bedrock bottom layers, are derived from a number of datasets, namely the Food and Agriculture Organization soil database (FAO, 1988), the pan-European River and Catchment Database (Vogt et al., 2007), International Hydrogeological map of Europe (IHME) (Duscher et al., 2015), and GLobal HYdrogeology MaPS (GLHYMPS) (Gleeson et al., 2014). The information on subsurface aquifers is derived from IHME. The bedrock geology is constructed from the IHME hydrogeological information and the lower-resolution GLHYMPS, in combination with the pan-European River and Catchment Database. The pan-European River and Catchment Database serves as a proxy for the alluvial aquifer system in ParFlow, which is assumed to lie underneath or near existing rivers.

2.3 Multi-model GCM–RCM ensemble

RCMs driven by CMIP5 GCMs control simulations (r1i1p1 ensemble members) over the European domain from the CORDEX experiment at EUR-11 horizontal resolution are used in conjunction with the coupled TSMP modelling platform to study the impact of 3D groundwater dynamics on the statistics of heat events. Note that CMIP5 GCM historical control simulations are performed under observed natural and anthropogenic forcing (e.g., Taylor et al., 2012).

Giorgetta et al. (2013)Voldoire et al. (2013)Bentsen et al. (2013)Dunne et al. (2012)Dufresne et al. (2013)Shrestha et al. (2014)Rockel et al. (2008)Pothapakula et al. (2020)Jacob and Podzun (1997)

Table 1 The GCM–RCM climate change scenario control runs.

Based on data availability, this study considers the following CORDEX ensemble members identified by their institutions: CLMcom (CCLM4-8-17 driven by MPI-ESM-LR and CNRM-CM5), CLMcom-ETH (COSMO-crCLIM driven by MPI-ESM-LR, CNRM-CM5, and NCC-NorESM1-M), MPI-CSC (REMO2009 driven by MPI-ESM-LR), and GERICS (REMO2015 driven by NCC-NorESM1-M, NOAA-GFDL-ESM2G, and IPSL-CM5A-LR); see Table 1 for details. Such a multi-model GCM–RCM ensemble includes two main groups of RCMs, namely COSMO and REMO in different versions, and five different GCMs, for a total of 10 different GCM–RCM pairs. TSMP is most compatible with CCLM4-8-17, with the main differences in the COSMO lower boundary condition: in TSMP, the lower boundary condition accounts for groundwater feedbacks due to the coupling between the CLM land surface model and the ParFlow hydrological model, unlike in CCLM4-8-17, where soil processes are modelled by TERRA-ML, the soil–vegetation land surface model of COSMO (e.g., Grasselt et al., 2008; Schlemmer et al., 2018). With the exception of TSMP, the RCMs in the considered ensemble lack the closure of water and energy cycles due to simplifications of the representation of subsurface hydrodynamics.

2.4 Detection and analysis of heat events

There is no universally accepted method for detecting heat events, but the most commonly used approach is based on a percentile temperature threshold (e.g., Zhang et al., 2005, 2011; Sulikowska and Wypych, 2020). However, temperature-based diagnostics are often ambiguous or inconsistent and only partially describe heat events (Perkins and Alexander, 2013).

Figure 1 Schematic of a summer heat wave (HW) detection. An example is given for June–July–August of 1972 for one grid element [250, 300] of the European CORDEX domain. Data taken from the TSMP simulation. The solid black line is the daily mean 2 m air temperature for the summer season of 1972. The dashed green line is the climatological daily mean 2 m air temperature for 1961–1990, and the solid green line shows its smoothing using the Butterworth filter. The solid violet line represents the 90th percentile of the daily mean 2 m air temperature calculated from a 5 d window centred on each summer calendar day of the 1961–1990 reference period. The shaded light red colour indicates days with temperatures above the climatological mean, and the shaded dark red colour highlights days with temperatures above the 90th percentile. The characteristics of the heat events (start and end dates, duration, intensity) detected during the considered summer season are also given.

In this study, we determine a day with a daily mean temperature above the local 90th percentile of the reference period as a hot day. We calculate the 90th percentile for every summer day and for each EUR-11 grid point of the European CORDEX domain from a consecutive 5 d moving window centred on that calendar day from the 30-year period between 1961 and 1990 in each RCM considered. The first occurrence of a hot day defines the beginning of a heat event. A series of hot days constitutes a heat event, highlighted in dark red in Fig. 1. A heat event is interrupted if the daily mean temperature drops below the 90th-percentile-based threshold.

A heat event can be characterized by its duration, intensity, and frequency (e.g., Horton et al., 2016). A heat event duration is the number of consecutive days over which the heat event lasts. If a heat event lasts long enough, it can be classified as a heat wave. Similar to Fischer and Schär (2010), we define a heat wave as a spell of at least 6 consecutive days with daily mean 2 m air temperatures above the local 90th percentile of the 1961–1990 reference period; see Fig. 1. The total number of hot days during the investigated period corresponds to the TG90p heat index from the joint CCl/CLIVAR/JCOMM Expert Team on Climate Change Detection and Indices (ETCCDI) (Alexander et al., 2006) and describes the number of days with TG_ij>TG_in90, where TG_ij is a daily mean temperature on day i of the investigated period j, and TG_in90 is the 90th percentile calculated for day i from the 30-year reference period n. Note that we consistently use the terms hot day, heat event, and heat wave throughout this article.

A heat event intensity is the maximum of the difference between the daily mean temperature and the 90th percentile of the reference period within a single heat event (e.g., Vautard et al., 2013). Intensity represents the severity of a heat event (see Fig. 1). Adopting the definition of the heat wave duration index (HWDI) from Frich et al. (2002), here we classify a heat wave as intense if its intensity is at least 5 K. In literature, there are other classifications of heat waves depending on their intensity, e.g., low intensity, severe, and extreme (Nairn and Fawcett, 2014) or weak, moderate, and intense (Lhotka and Kyselý, 2015).

The frequency of heat events of a certain type (e.g., of a certain duration or intensity) over the investigated period is the number of these heat events divided by the total number of all heat events that occurred during the period under study (e.g., Vautard et al., 2013). For example, in Fig. 1, heat events with a duration of 2 d occur four times, and the total number of all heat events is equal to 10. Therefore, the resulting frequency of 2 d heat events during that summer season is 0.4, indicating that 40 % of all heat events are 2 d in duration.

3 Results

In the following we examine the impact of groundwater representation on the distribution of simulated heat events in RCMs over Europe. In particular, we investigate whether the new dataset from TSMP driven by MPI-ESM-LR is consistent with the GCM–RCM climate change scenario control runs from the CORDEX ensemble (see Table 1), and we seek to gain insight into the role of explicit groundwater representation in long-term climate simulations. We assess the statistics of the characteristics of heat events, i.e., their duration, intensity, and frequency, from daily mean 2 m air temperatures on the native EUR-11 grid. The analysis is conducted in the focus domain (Fig. 2), which covers the European continent (36–70° N, 10° W–30° E). The analysis is performed for the summer seasons of the 30-year period from 1976–2005 relative to the reference period 1961–1990 for each RCM. Note that grid elements belonging to the ocean are omitted from the analysis.

3.1 Number of hot days

We examine the impact of groundwater dynamics on the interannual variability of the occurrence of hot days during the summer seasons from 1976–2005 regarding the reference period 1961–1990 in each RCM. A comparison of the time series of the number of hot days, i.e., the TG90p index, averaged over the summer season and over the focus domain shows that the impact of an explicit representation of groundwater dynamics in RCMs varies from year to year (Fig. 3). The summer TG90p index averaged between 1976 and 2005 results in 10.95 d in the TSMP simulation and 11.64 d in the CORDEX multi-model ensemble mean (see Fig. 3). A positive linear trend in the summer TG90p index is observed in all considered RCMs, with the decadal change in the TSMP simulation being 1.53 d, while this value averaged over the CORDEX ensemble reaches 1.99 d.

The spatial distributions of the seasonal mean, variability, and decadal change in the summer TG90p index are shown in Figs. 4–6. There, the spatial patterns from RCMs driven by the same GCMs show rather similar behaviour, indicating that the climatological occurrence of summer hot days is largely controlled by the large-scale atmospheric circulation imposed by the GCM boundary conditions. TSMP produces the smoothest spatial distribution of the seasonal mean and variability of the summer TG90p index compared to the CORDEX ensemble (see standard deviations in Figs. 4, 5). The mean and interannual variability of the summer TG90p index averaged over the focus domain are also lowest in TSMP compared to the CORDEX ensemble (Tables A1 and A2 in Appendix A).

The TSMP-simulated summer TG90p index is consistent with that of the RCMs driven by MPI-ESM-LR from the CORDEX ensemble, although there are some regional differences (see Figs. 4a–d, 5a–d, 6a–d). The largest differences occur in the Iberian Peninsula PRUDENCE region, with TSMP yielding the lowest values: the mean TG90p index is equal to 10.36 d in TSMP and ranges from 12.54 to 12.75 d in the CORDEX RCMs driven by MPI-ESM-LR, variability of the TG90p index reaches 6.17 d in TSMP and ranges from 8.01 to 9.59 d in the CORDEX RCMs driven by MPI-ESM-LR, and decadal change in the TG90p index is 2.26 d in TSMP and ranges from 3.66 to 4.25 d in the CORDEX RCMs driven by MPI-ESM-LR (see Appendix A). As for the decadal change, the RCMs driven by MPI-ESM-LR show a positive trend in southern and central Europe and a negative trend in northern Europe. Note that there is no unequivocal agreement in the decadal change at the entire GCM–RCM ensemble considered, yet the decadal trend of the summer TG90p index averaged over the focus domain is positive in all models of the investigated GCM–RCM ensemble, reaching values between 1.13 and 2.71 d (see Table A3 in Appendix A).

Figure 2 Mean 90th percentile of 2 m air temperatures for the summer season, based on TSMP simulation from 1961 to 1990. The white box indicates a focus domain (36–70° N, 10° W–30° E) used in the analysis. PRUDENCE regions are shown with grey boxes: British Isles (BI), Iberian Peninsula (IP), France (FR), mid‐Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA).

Figure 3 Time series and linear trends of the summer TG90p index, averaged over the focus domain, during 1976–2005 in relation to the reference period 1961–1990, in the TSMP simulation and the CORDEX ensemble. The solid and dashed red lines show the summer mean TG90p index and its linear trend from the TSMP simulation. The black and grey lines represent the summer mean TG90p index from the GCM–RCM of the CORDEX ensemble and the green lines are their linear trends, respectively. The summer TG90p index averaged over the CORDEX multi-model ensemble is shown with the solid blue line, and its linear trend is shown with the dashed blue line.

3.2 Heat events of different durations

The summer seasonal number of heat events (i.e., series of consecutive hot days) of different durations occurring on average between 1976 and 2005 in the focus domain in TSMP and the CORDEX ensemble is shown in Fig. 7a. In the considered GCM–RCM multi-model ensemble, the mean number of heat events of any duration per summer in the focus domain ranges from 4.18 to 4.86 heat events, with TSMP simulating 4.66 heat events. The ratios of the number of heat events between the CORDEX RCMs and TSMP (see blue lines in Fig. 7a) are greater than 1 for heat waves, i.e., heat events lasting at least 6 d, and increases towards the heat waves of long durations. This behaviour indicates that TSMP systematically simulates the least number of heat waves compared to the CORDEX ensemble. An intercomparison of the RCMs from the CORDEX ensemble shows that COSMO tends to simulate fewer heat waves on average over the focus domain compared to REMO.

The GCM–RCM multi-model ensemble leads to different spatial distributions of heat waves (Fig. 8). TSMP generates the smoothest spatial distribution of the number of heat waves, resulting in the smallest regional differences compared to the CORDEX ensemble (see indicated standard deviations in Fig. 8). The decadal number of summer heat waves on average over the focus domain ranges from 3.25 to 5.09 in the GCM–RCM multi-model ensemble, with the lowest value in TSMP (Table B3 in Appendix B). When comparing TSMP with the CORDEX RCMs driven by MPI-ESM-LR, the TSMP simulation has most of the heat waves located in central Europe, while the RCMs from the CORDEX ensemble simulate the highest number of heat waves towards southern Europe. The largest differences occur in the Iberian Peninsula, with 2.51 heat waves per decade in TSMP and between 4.88 and 5.25 heat waves per decade in the RCMs driven by MPI-ESM-LR from the CORDEX ensemble.

Figure 4 Spatial distribution of the summer TG90p index averaged between 1976 and 2005 for TSMP (a) and the CORDEX ensemble (b–j). A standard deviation (SD) of the spatial distribution of the summer TG90p is indicated in every panel.

Figure 5 Variability of the summer TG90p index, calculated for each grid element from the time series of the summer mean TG90p index between 1976 and 2005 as standard deviation, for TSMP (a) and the CORDEX ensemble (b–j). A standard deviation (SD) of the spatial distribution of the variability of the summer TG90p index is indicated in every panel.

Figure 6 Spatial distribution of the decadal change in the summer TG90p index, calculated for each grid element from the time series of the summer mean TG90p index between 1976 and 2005 as a linear trend, for TSMP (a) and the CORDEX ensemble (b–j).

The contribution of heat waves to the total number of hot days is presented in Fig. 9. Here, heat waves account for from 22.38  % to 34.40  % of hot days on average in the focus domain, with TSMP giving the lowest value (Table B2 in Appendix B). The highest fraction of hot days attributed to heat waves prevails in Scandinavia, from 26.50  % to 39.04  %, while eastern Europe tends to be the region with the lowest number of hot days associated with heat waves. From the comparison of TSMP with the RCMs driven by MPI-ESM-LR, the largest differences in the proportion of hot days belonging to heat waves are observed in the Iberian Peninsula, where TSMP simulates 17.47  % and the RCMs driven by MPI-ESM-LR from the CORDEX ensemble simulate between 31.19  % and 33.04  %.

Figure 7 (a) Mean number of summer heat events (HEN, y axis) of duration equal to or greater than a given number of days (x axis) as a function of this number of days. The averaging is performed over the focus domain and the total number of investigated years, i.e., 30 years, from 1976–2005. HEN is shown with the solid red line for TSMP and with the black and grey lines for the CORDEX ensemble. The total HEN occurring on average over the summer season in the focus domain in each RCM considered is given in the inset table. (b) Frequency of heat waves (HWF, y axis) with intensities equal to or higher than a value indicated on the abscissa, occurring in the focus domain from 1976–2005, as a function of this intensity. HWF is shown with the red solid line for TSMP and with the black and grey lines for the CORDEX ensemble. Panels (a) and (b) also show the ratios between the RCMs from the CORDEX ensemble and TSMP for HEN and HWN values, and are represented by the blue lines, respectively. Data are taken from the summer seasons between 1976 and 2005 relative to the reference period 1961–1990 in each RCM. The representation of the dependencies is adopted from the work of Vautard et al. (2013).

Figure 8 Spatial distribution of the decadal number of summer heat waves (HWN), calculated from the data between 1976 and 2005 relative to the reference period 1961–1990, for TSMP (a) and the CORDEX ensemble (b–j). A standard deviation (SD) of the spatial distribution of the decadal HWN is indicated in every panel.

Figure 9 Contribution of heat waves to the number of hot days (%), calculated from the total number of heat waves and hot days accumulated between 1976 and 2005, for TSMP (a) and the CORDEX ensemble (b–j).

Figure 10 Frequency of intense heat waves (HWF), i.e., heat waves with intensities of at least 5 K, relative to the total number of heat waves occurring between 1976 and 2005 in each RCM. The HWF distribution is shown for TSMP (a) and the CORDEX ensemble (b–j).

3.3 Heat waves of different intensities

The dependencies of the frequency of heat waves that occurred in the focus domain in the period from 1976–2005 on their intensities are shown in Fig. 7b. The maximum frequency of heat waves is equal to 1 for an intensity greater than 0 because all heat waves are taken into account for each RCM. The ratios of heat wave frequencies between the CORDEX RCMs and TSMP (see blue lines in Fig. 7b) increase towards intense heat waves, i.e., heat waves with an intensity of at least 5 K, except for REMO2015 driven by IPSL-CM5A-LR. TSMP shows a systematic behaviour to simulate less intense heat waves on average in the focus domain compared to the CORDEX ensemble. Note that the largest discrepancy in the frequency of heat waves with TSMP is found in CCLM4-8-17 driven by MPI-ESM-LR, up to a factor of 12 or even more depending on the intensity considered, although TSMP is the most compatible with CCLM4-8-17 in the considered GCM–RCM multi-model ensemble. An intercomparison of the RCMs within the CORDEX ensemble shows that COSMO tends to simulate more intense heat waves than REMO on average over the focus domain.

The most intense heat waves are located in western and northern Europe in the majority of RCMs of the considered multi-model ensemble (Fig. 10). The frequency of intense heat waves occurring in the focus domain between 1976 and 2005 ranges from 0.174 to 0.301, i.e., from 17.4 % to 30.1 % heat waves are intense in the GCM–RCM multi-model ensemble, with TSMP giving the second lowest value after REMO2015 driven by IPSL-CM5A-LR (Table B3 in Appendix B). As already noted, TSMP has the largest discrepancy in the frequency of heat waves with CCLM4-8-17 driven by MPI-ESM-LR, with particularly large differences in the France PRUDENCE region, where TSMP leads to a frequency of 0.246 and CCLM4-8-17 to 0.468. It is important to point out that the regions with the highest number of intense heat waves do not necessarily coincide with the regions that experience the most heat waves. The origin of such behaviour should be further investigated and is beyond the scope of this analysis.

4 Discussion

4.1 Physical mechanisms

Compared to the simplified 1D free-drainage approach, the physics-based 3D groundwater representation in TSMP results in regionally shallower groundwater levels, causing wetter soils (Keune et al., 2016). This leads to enhanced evapotranspiration by an increase in the latent heat flux and a decrease in the sensible heat flux (e.g., Maxwell and Condon, 2016). In turn, higher evapotranspiration causes moistening of the lower atmosphere and increases downward longwave radiation due to the greenhouse effect of water vapour; on the other hand, it causes cooling of the surface and reduces outgoing surface longwave radiation (e.g., Pal and Eltahir, 2001; Yang et al., 2018). Additionally, higher evapotranspiration may lead to moist convection or rainfall and further affect soil moisture. The simplified representation of groundwater dynamics in RCMs leads to the opposite effect, i.e., deeper groundwater levels and drier soils, overestimating the coupling between the land surface and the atmosphere. This causes a decrease in cloud cover and an enhancement of net solar radiation, thus increasing near-surface temperatures, which further reduces soil moisture (e.g., Vogel et al., 2018; Hartick et al., 2022).

The results of our study suggest that the response of summer heat events to the physics-based 3D representation of groundwater in TSMP varies from year to year (see Fig. 3). Unlike the RCMs from the CORDEX ensemble, the explicit representation of groundwater incorporated into the TSMP model accounts for long-term soil moisture memory effects (Hartick et al., 2021). Soil moisture memory contributes to either increasing the probability of a subsurface water storage deficit in regions that have had a subsurface water deficit in the previous year due to drought conditions, thereby increasing the occurrence of heat events, or, conversely, buffering droughts and reducing the number of heat events (e.g., Martínez-de la Torre and Miguez-Macho, 2019; Hartick et al., 2021; Dirmeyer et al., 2021). Droughts can also remotely affect areas outside the drought region through changes in atmospheric circulation and advection of air masses and further contribute to the evolution of heat events (Fischer et al., 2007).

Considering an extended period of 30 years, from 1976–2005, TSMP driven by MPI-ESM-LR shows systematic differences in the distribution of heat events characteristics (i.e., duration, intensity, frequency) compared to the CORDEX ensemble by simulating fewer, shorter, and less severe heat events in Europe with smaller regional differences. We relate this behaviour to a more realistically simulated soil moisture and thus evapotranspiration in TSMP (see also Furusho-Percot et al., 2022). In addition, the tendency for different responses in the different PRUDENCE regions can be explained by the soil moisture–temperature feedbacks associated with evaporative regimes: energy-limited in northern Europe (i.e., a wet regime in which land evaporation is mainly controlled by incoming radiation) and moisture-limited prevailing in southern Europe (i.e., a dry regime in which land evaporation increases or decreases in response to increased or decreased soil moisture content) (e.g., Seneviratne et al., 2010; Haghighi et al., 2018; Jach et al., 2022). The comparison of TSMP with the most compatible RCM from the CORDEX ensemble, i.e., CCLM4-8-17 forced by MPI-ESM-LR, shows that groundwater representation has a particularly strong influence on the intensity of heat waves versus their duration (see Fig. 7); the physical mechanisms of this phenomenon require further investigation and are beyond the scope of this study.

4.2 Methodology limitations

In order to capture the full range of divergence in simulation results over the historical period, and hence potential uncertainties in the climate projections, it is necessary to combine as many different GCM–RCM pairs as possible (e.g., Déqué et al., 2012; Christensen and Kjellström, 2020). Often, certain RCMs and GCMs are overrepresented compared to others, even leading to contradictions in climate change signals (Turco et al., 2013; Fernández et al., 2019). Note that the role of the GCM-imposed boundary conditions is generally greater that the role of RCMs (Déqué et al., 2007; Evin et al., 2021). In this study, a limited number of GCM–RCM pairs (see Table 1) were used for comparison with TSMP; therefore, expanding the multi-model GCM–RCM ensemble or using other GCM–RCM pairs may lead to quantitatively different results.

It is important to emphasize that the GCM–RCM multi-model ensemble is not intended for direct comparison between individual ensemble members, as it includes different RCMs in combination with different driving GCMs. Due to the interplay of various factors (e.g., conceptual and structural differences between models, model set-up, different physical parameterizations, internal variability, boundary conditions, representation of the subsurface–land–atmosphere feedbacks) in addition to groundwater representation, it is challenging to reveal the exact cause–effect relationships between the explicit groundwater representation and simulated heat events. However, considering an extended period of 30 years allows us to draw statistical conclusions.

To quantify the exact impact of the explicit representation of groundwater in TSMP and minimize the impact of other factors, it would be necessary to additionally conduct a long-term TSMP climate simulation with a simplified 1D free-drainage approach for groundwater representation and then compare the affected processes within TSMP rather than across the multi-model GCM–RCM ensemble. Since our study uses the same version of the TSMP model as in Keune et al. (2016), which have already shown the effect of 3D groundwater dynamics on the water and energy balance, and taking into account the high computational cost, the additional TSMP simulation with simplified groundwater representation is not feasible. The multi-model GCM–RCM ensemble analysis performed in this study provides valuable insights into the consistency between the new dataset of TSMP driven by MPI-ESM-LR and the CORDEX ensemble.

Note that the results of this research are limited to the definitions used (e.g., hot day, heat event, heat wave, intense heat wave), the method of percentile estimation, and the choice of the investigated and reference time periods. For instance, Sulikowska and Wypych (2020) indicate that different variations in the metrics lead to different distributions of hot days during the summer season in Europe.

5 Summary and conclusions

We presented a first-of-its-kind dataset of TSMP simulation driven by the CMIP5 MPI-ESM-LR GCM boundary conditions, in the context of dynamical downscaling of GCMs by RCMs for climate change studies. Unlike most RCMs, TSMP is a fully coupled regional climate system model with an explicit representation of groundwater, enabling complete water and energy cycles to be considered. We examined the role of groundwater representation for heat events in a multi-model ensemble of 10 different GCM–RCM pairs by comparing TSMP results and those from the CORDEX experiment with simplified groundwater representation. Specifically, we performed a statistical analysis of the characteristics of heat events (i.e., duration, intensity, frequency) over Europe during the summer seasons between 1976 and 2005 in relation to the reference period of 1961–1990.

The characteristics of heat events simulated by TSMP are consistent with the CORDEX ensemble, although systematic differences are observed, which we attribute to an explicit representation of groundwater in TSMP. Our findings suggest that the integrated 3D physics-based representation of groundwater in TSMP leads to a reduction in the number of hot summer days, their interannual variability, and decadal change and causes smaller regional differences compared to the CORDEX ensemble. The explicit representation of groundwater also influences simulated heat waves, resulting in a decrease in the number of heat waves and a reduction in their duration and intensity compared to the CORDEX ensemble. From the comparison of TSMP and the CORDEX RCMs driven by MPI-ESM-LR, the Iberian Peninsula is the most sensitive region with respect to groundwater representation.

This study clearly indicates that a coupled regional climate system model with 3D groundwater dynamics systematically simulates a different climatology of heat events in Europe compared to RCMs with simplified representation of groundwater. The results emphasize the importance of hydrological processes for reliable climate modelling and for reducing biases in the duration, intensity, and frequency of simulated heat waves in particular, and they are of further importance when assessing uncertainties in climate change projections.

Appendix A: TG90p index in different regions of Europe

Table A1 Summer TG90p index (d) averaged between 1976 and 2005 in the GCM–RCM multi-model ensemble for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 4 for the spatial distribution.

Table A2 Variability of the summer TG90p index (d), calculated from the data between 1976 and 2005 in the GCM–RCM multi-model ensemble, for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 5 for the spatial distribution.

Table A3 Decadal change in the summer TG90p index (d), calculated from the data between 1976 and 2005 in the GCM–RCM multi-model ensemble, for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 6 for the spatial distribution.

Appendix B: Characteristics of heat waves for different regions of Europe

Table B1 Decadal number of summer heat waves, calculated from the data between 1976 and 2005 in the GCM–RCM multi-model ensemble, for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 8 for the spatial distribution.

Table B2 Contribution of heat waves to the number of hot days (%), based on the data from 1976–2005 in the GCM–RCM multi-model ensemble, for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 9 for the spatial distribution.

Table B3 Frequency of intense heat waves, calculated from the data between 1976 and 2005 in the GCM–RCM multi-model ensemble, for the focus domain (FD; see Fig. 2) and the PRUDENCE regions: British Isles (BI), Iberian Peninsula (IP), France (FR), mid-Europe (ME), Scandinavia (SC), Alps (AL), Mediterranean (MD), and eastern Europe (EA). Refer to Fig. 10 for the spatial distribution.