2.1 Scope and scale of the study

The scope of the study covers mainland France from January 1, 2010, to December 31, 2019. The study of the relationship between mortality and environmental variables faces an inherent trade-off between bias and variance. This trade-off reflects the discretization required by the analysis. On the one hand, the finer the spatial and geographical scale of the study, the more the environmental data reflect the reality locally experienced. Spatial or temporal averages smooth out extreme values. On the other hand, if the scale is too small, the number of deaths observed becomes too small, and the statistical robustness is lost. A 5-day interval was chosen, at a small regional scale between cities and departments. The construction of these areas is described below in Section 2.2. This spatial grid differs from those used in the two most closely related studies. In Robben et al. (Reference Robben, Antonio and Kleinow2025), the spatial resolution is departments (corresponding to the Nomenclature of Territorial Units for Statistics (NUTS) of level 3), and in Guibert et al. (Reference Guibert, Pincemin and Planchet2025), temperatures are averaged for all mainland France. This territorial division seems to be too large. Indeed, choosing to average temperatures over France eliminates the extremes that can occur very locally, for example, on August 27, 2015, 34 $^\circ$ C was measured in Dax and 14 $^\circ$ C in Chartres. The average temperature was 24 $^\circ$ C, while the south-west experienced scorching temperatures. Moreover, both territorial divisions make it difficult to take into account the UHI effect. Due to the modification of land surfaces and the waste heat generated by energy usage, temperatures in urban areas can be much higher than in surrounding rural areas, particularly at night. This effect has a significant impact on mortality (Wong et al. (Reference Wong, Paddon and Jimenez2013), Yadav et al. (Reference Yadav, Rajendra, Awasthi, Singh and Bhushan2023), among others) and must therefore be taken into account. The time scale also affects the results and must therefore be chosen in conjunction with the spatial scale. An appropriate temporal scale must be sufficiently fine, applicable in the same way whatever the year, and take into account the effects on mortality associated with the day of the week. Three scales appear to be particularly relevant: daily, as in Guibert et al. (Reference Guibert, Pincemin and Planchet2025), weekly, as in Robben et al. (Reference Robben, Antonio and Kleinow2025), and the 5-day scale. Given the fine spatial grid adopted, a daily time step was not feasible. A particular feature of mortality is its dependence on the day of the week. A slightly higher proportion of deaths is reported on Mondays and Tuesdays. Yet the effect is very small. Indeed, in the study data, 14.5% of deaths occur on Tuesdays and 13.9% on Sundays. When considering rolling 5-day periods, the maximum difference between the maximum number of deaths in a period and the minimum number is 0.8%: 71.9% of deaths occurred from Monday to Friday, compared with 71.2% from Wednesday to Sunday. Therefore, considering 5-day periods introduces a negligible bias. The advantage is that 5 divides a 365-day year into equal parts. February 29 in leap years has not been considered. This event only occurs twice out of the 3,652 days considered in the study, which introduces a negligible bias. Deaths are treated according to sex and age. Ages are grouped in 5-year increments, from 40 to 99.

2.2 Data

Deaths data. The mortality data come from the official death register maintained by INSEE. Since 1970, INSEE has compiled individual-level information on all deaths reported by French municipalities as part of their public service obligations. For each deceased individual, the dataset records sex, exact dates of birth and death, and the place of birth and death, identified by the official commune codes. The level of detail is therefore highly granular: ages are available on a continuous scale (not in intervals), places can be traced to the commune-level, and no aggregation is applied at the source. Monthly files are released and subsequently concatenated into annual and decennial compilations, which are the “ten-year files” used in this study. These files are not sample-based but exhaustive and cover the entire French population, including deaths occurring abroad among French nationals. For our analysis, we restricted the scope to deaths occurring in mainland France, excluding overseas territories.

Table 1.

Climatic data selected from Météo France for the study

Figure F.1 Long description

The scatter plot compares observed and predicted mortality ratios for training and testing sets. The x-axis represents the observed over baseline realization, and the y-axis represents the predicted over baseline expectation. Small semi-transparent points represent individual 5-day period realizations, reflecting high variance in death counts. Large bordered circles represent binned means for 50 quantiles, showing the model’s expected value against the average realization. Color-coding indicates mean temperature in degrees Celsius. The dashed diagonal line represents ideal calibration where the model expectation perfectly matches the observed mean. The plot shows a comparison between observed and predicted mortality ratios, with clusters and patterns indicating the model’s performance. All values are approximated.

Exposure data. Population exposures are obtained from official INSEE files that report resident population by single year of age, sex, and commune (excluding Mayotte) from 2006 onwards, with additional subdivisions for the districts of Paris, Lyon, and Marseille. These files provide annual cross-sections of the population as of January 1 of each year. To construct exposures at finer temporal resolutions, we follow standard actuarial practice and assume that the population evolves linearly between successive annual counts. The calculation is described by Equation 3.

Climate data. Since the beginning of 2024, Météo France has provided online access to a vast amount of weather and climate data. Specifically, it provides comprehensive climatological data from all the stations in metropolitan France and overseas from the beginning of the recording period for all the available parameters. Climatological checks have been carried out on these data. Two types of data are available for this type of study: on the one hand, in situ data, corresponding to measured and observed data. On the other hand, reanalysis data obtained by combining past observations with models to generate consistent time series of multiple climate variables. In the case of a small number of observations, reanalysis data are highly advantageous, especially as it has been shown that reanalysis data can give good results in constructing a relationship between temperature and mortality (de Schrijver et al., Reference de Schrijver, Folly, Schneider, Royé, Franco, Gasparrini and Vicedo‐Cabrera2021). In the present case, a large and sufficient amount of in situ high-quality data was available as illustrated in Figure A.1 in Appendix A, which shows the very fine spatial grid of the weather stations with open-source data selected for this study. Of all existing weather stations in the data, 501 were finally selected, based on the percentage of missing values in the main weather variables. In total, the Météo France data consists of 66 fields. However, many of these fields are incomplete. We have selected a limited number of climate variables for each weather station. Data are described in Table 1. Of these stations, only 0.65% of values are missing. The variable with the most missing data is the wind speed, with an occurrence of 1.50%. Missing values are replaced by linear interpolations of data from the nearest days for the same station.

Meteorological data were linked to other environmental covariates, specifically those pertaining to air pollution and vegetation. Adélaide et al. (Reference Adélaïde, Hough, Seyve, Kloog, Fifre, Launoy, Launay, Pascal and Lepeule2024) previously analyzed the correlation between heat exposure, air pollution, lack of vegetation, and social deprivation from a public health perspective. Accordingly, pollution and vegetation are regarded here as contributing factors that interact with climatic stressors to influence mortality outcomes. Consequently, our study incorporates data on pollution and vegetation.

Pollution data. Air pollution is a well-established determinant of mortality. Pascal et al. (Reference Pascal, Falq, Wagner, Chatignoux, Corso, Blanchard, Host, Pascal and Larrieu2014) demonstrated the relationship between mortality and particulate matter (PM) in France. Furthermore, high ozone concentrations are positively correlated with excess mortality in both the short (Vicedo-Cabrera et al., Reference Vicedo-Cabrera, Sera, Liu, Armstrong, Milojevic, Guo, Tong, Lavigne, Kyselý, Urban, Orru, Indermitte, Pascal, Huber, Schneider, Katsouyanni, Samoli, Stafoggia, Scortichini and Gasparrini2020) and long term (Atkinson et al., Reference Atkinson, Butland, Dimitroulopoulou, Heal, Stedman, Carslaw, Jarvis, Heaviside, Vardoulakis, Walton and Anderson2016). The data for this study were sourced from the French Institut national de l’environnement industriel et des risques (Ineris). These data are available at municipal and annual resolutions and were subsequently averaged to the regional and temporal scales of this study. The dataset comprises several variables, with those selected for analysis described in Table 2.

Land use and vegetation data. The CORINE Land Cover (CLC) is a biophysical inventory of land cover, produced through the visual interpretation of satellite images in accordance with a specific nomenclature. Data can be accessed through the website of the French Ministry of Environmental Transition, which collects statistics and data pertaining to sustainable development. The 15-level land cover dataset for mainland France, produced in 2018, was selected. The fifteen are the following: urban fabric; industrial, commercial, and transport units; mine, dump, and construction sites; artificial, nonagricultural vegetated areas; arable land; permanent crops; pastures; heterogeneous agricultural areas; forests; scrub and/or herbaceous vegetation associations; open spaces with little or no vegetation; inland wetlands; coastal wetlands; inland waters; marine waters. The data are available at the municipality level. For each municipality, the surface area of each of these categories is available. Data on the geographical scale of the study were aggregated and then normalized so that the sum for a geographical area is 1.

Climate type data. France’s distinctive geographical location and topography result in a diverse range of climatic conditions within a relatively confined geographical area. France is bordered by two seas and an ocean and is characterized by the presence of several mountainous areas. Joly et al. (Reference Joly, Brossard, Cardot, Cavailhes, Hilal and Wavresky2010) proposed a typology of climate together with a statistical method for its construction. A multivariate analysis of the climatic variables was employed. The resulting classification identifies a subdivision of France into eight types.

Construction of geographical scale of the study. The available data permit mortality studies at various geographical scales, ranging from the commune-level (34,935 communes as of January 1, 2024) to the national level. We selected an intermediate scale: the catchment area (bassin de vie), as defined by INSEE. These areas represent the smallest territories in which residents have access to the most common facilities and services. To ensure compatibility with our study period, we utilized the 2012 definition, which comprises 1,666 catchment areas. As many of these areas are too sparsely populated for robust mortality analysis, we constructed larger zones that are both environmentally homogeneous and statistically significant.

Table 2.

Pollution data from Ineris selected for the study

Table 2. Long description

The table presents pollution data variables and their descriptions used in a study on air pollution in France. It includes six rows and two columns. The first column lists variable symbols such as Mean_NO2, Mean_O3, Mean_somo35, Mean_AOT40, Mean_PM10, and Mean_PM2_5. The second column provides descriptions of these variables, detailing what each represents, such as average concentrations of NO2, O3, PM10, and PM2.5 over specific periods and conditions. The table helps in understanding the different pollution metrics analyzed in the study.

To maintain geographical coherence, we constructed path-connected clusters using agglomerative hierarchical clustering. In this approach, each observation begins in its own cluster, and pairs are progressively merged. We employed the Agglomerative Clustering implementation from the Python library scikit-learn. Each catchment area is described by a feature vector $x \in \mathbb{R}^{34}$ comprising: population size (1 variable), geographical descriptors (latitude, longitude, altitude: 3 variables), pollution indicators (6 variables, see Table 2), land use shares (15 variables), and climate type dummies (9 variables). All variables were standardized prior to clustering to ensure each dimension contributed equally to the distance measure. We applied Ward’s minimum-variance criterion with Euclidean distance while enforcing spatial contiguity through a Queen contiguity matrix $W$ , which restricts merges to adjacent catchment areas. Ward’s method minimizes the total within-cluster variance; at each step, the algorithm merges the pair of clusters that results in the minimum increase in total within-cluster variance (WCV). Formally, the algorithm merges clusters A and B that minimize:

(1) \begin{equation} \min _{A,B} \sum _{x\in A\cup B} \lVert x - \mu _{A\cup B} \rVert ^2 - \sum _{x\in A} \lVert x - \mu _{A} \rVert ^2 - \sum _{x\in B} \lVert x - \mu _{B} \rVert ^2 \end{equation}

where $\mu _X$ denotes the centroid of cluster $X$ . This procedure yielded an initial partition of 150 clusters.

A pruning step was subsequently applied to address small residual clusters with fewer than 100,000 inhabitants. Any such cluster $r$ was reassigned to the most similar neighboring cluster j based on Euclidean distance in the feature space:

(2) \begin{equation} j^* = \arg \min _{j \in \mathcal{N}(r)} \lVert \mu _r - \mu _j \rVert , \end{equation}

where $\mathcal{N}(r)$ is the set of adjacent clusters. The final partition consists of 89 geographically coherent and environmentally homogeneous zones. These are more suitable for mortality analysis than standard administrative divisions as they effectively capture key ecological and climatic factors. Large-area clusters typically correspond to rural zones, whereas small-area clusters represent dense urban environments. Figure 1 shows the clustering methodology: Step 1 shows the 150 clusters obtained after the initial algorithm; Step 2 identifies the “residual” clusters (those below the population threshold) designated for pruning; Step 3 displays the final 89 clusters used in the analysis.

Figure 1

Step 1 shows the 150 clusters of mainland France. Step 2 illustrates the smallest clusters with a number of inhabitants lower than 100.000. Step 3 represents the final 89 clusters of mainland France.

Figure 1 Long description

The image consists of three maps of mainland France. Step 1 shows the initial division of France into 150 clusters. Step 2 highlights the smallest clusters, each with fewer than 100,000 inhabitants. Step 3 presents the final division, resulting in 89 clusters. Each step is visually distinct, with different color coding to represent the clusters.

2.3 Main notation

We adopt the notation used by Robben et al. (Reference Robben, Antonio and Kleinow2025). To simplify the notation, we omit an index for sex; however, the analysis was disaggregated by sex, and the subsequent results reflect this distinction. Let $d_{x,t,p}^{(r)}$ be the observed death count in region $r$ at age group $x$ during the 5-day period $p$ in year $t$ , $p \in \mathcal{P} = \{0, 1, \ldots 72\}$ , $t \in \mathcal{T} = \{2010, \ldots , 2019\}$ , $x \in \mathcal{X} = \{40, 45, 50, \ldots , 90, 95\}$ (we recall that we decided to consider age groups of 5 years), $r \in \mathcal{R} = \{1, 2, \ldots , R\}$ with $R=89$ the number or regions. Let $E_{x,t,p}^{(r)}$ denote the exposure, defined as the number of individuals alive, in region $r$ for age group $x$ during the 5-day period $p$ in year $t$ .

The exposure data described in Section 2.2 provide annual population counts $N_{x,t}^{(r)}$ for each age group and region. These counts were obtained by aggregating ages into five-year classes and municipalities into the regions defined via the clustering procedure described previously. Exposure for a given time step $p$ is calculated by assuming a linear evolution between successive annual counts:

(3) \begin{equation} E_{x,t,p}^{(r)} = \frac {(72 - p) \times N_{x,t}^{(r)} + p \times N_{x,t+1}^{(r)}}{72} \end{equation}

reflecting the 73 five-day periods that constitute a year, numbered from 0 to 72.

The crude central mortality rate, $m_{x,t,p}^{(r)}$ is defined as the ratio of observed deaths to the population exposure:

\begin{equation*} m_{x,t,p}^{(r)} = \dfrac {d_{x,t,p}^{(r)}}{E_{x,t,p}^{(r)}}. \end{equation*}

Following standard actuarial practice, we assume the force of mortality remains constant within each five-day interval. Under this piecewise constant assumption and the standard assumption that deaths follow a Poisson distribution (given the exposure), the Maximum Likelihood Estimate (MLE) of the force of mortality, $\hat {\mu }_{x,t,p}^{(r)}$ , is mathematically identical to the crude central mortality rate, $m_{x,t,p}^{(r)}$ (Pitacco et al., Reference Pitacco, Denuit, Haberman and Olivieri2009).

The probability that an individual in region $r$ and age group $x$ , alive at the start of period $p$ in year $t$ , will die within the subsequent five-day interval is denoted by $q_{x,t,p}^{(r)}$ , and its estimate $\hat {q}_{x,t,p}^{(r)}$ . Given the aforementioned assumption of a piecewise constant force of mortality, the estimate is derived from the functional relationship between the hazard and the survival probability:

\begin{equation*} {\hat {q}_{x,t,p}^{(r)}} \approx 1 - \exp \left (-\hat {\mu }_{x,t,p}^{(r)}\right )\!, \end{equation*}

for all $x\in \mathcal{X}$ , $t\in \mathcal{T}$ , $p \in \mathcal{P}$ and $r \in \mathcal{R}$ .

For any variable $x^{(r)}$ , x(r), the superscript $N$ (i.e. $x^{N}$ ) denotes national-level values. Further-more, where the period index $p$ is omitted, the notation refers to the annualized value.

3.1 Regional models for the annual baseline

We initially applied the standard Lee–Carter model (Lee & Carter, Reference Lee and Carter1992) at the national level. This model assumes that the logarithm of the central death rate $m_{x,t}$ for age $x$ in year $t$ is given by:

(4) \begin{equation} \ln\!(m_{x,t}) = \alpha _x + \beta _x \kappa _t + \epsilon _{x,t}, \end{equation}

where $\alpha _x$ represents the average age profile of mortality, $\kappa _t$ captures the time trend common to all ages, and $\beta _x$ reflects the sensitivity of mortality at age $x$ to changes in $\kappa ^{(r)}_t$ . The error term $\epsilon _{x,t}$ accounts for idiosyncratic variation. Parameters were estimated following the classical approach: first, we computed the logarithm of the central death rates, and subsequently applied singular value decomposition (SVD) to isolate the age and time components. We then evaluated the relative difference between observed deaths and those predicted by the national Lee–Carter mortality curve when applied uniformly to each region, as illustrated in Figure 2. The expected deaths for age group x and year t are calculated as $\widehat {d}^{\,(r)}_{x,t} \;=\; E^{(r)}_{x,t}\, \hat {m}_{x,t}$ . Figure 2 displays the relative deviation for each region $r$ :

(5) \begin{equation} \left (\sum _{x,t}\widehat {d}^{\,(r)}_{x,t} - \sum _{x,t}d^{(r)}_{x,t} \right )\quad \Big/ \quad \sum _{x,t}d^{(r)}_{x,t} \end{equation}

The results highlight that mortality levels vary substantially across the territory. Furthermore, the analysis identifies systematic discrepancies between rural and urban regions, confirming the necessity of incorporating regional heterogeneity when constructing baseline mortality models.

Figure 2

Relative difference between the number of deaths observed and the number of deaths predicted by the national mortality curve obtained using the Lee-Carter method.

Figure 2 Long description

A heat map of France displays the relative difference between the number of deaths observed and the number of deaths predicted by the national mortality curve using the Lee-Carter method. The map uses a color scale ranging from dark purple to yellow, where darker colors represent lower deviations and lighter colors indicate higher deviations. The map shows that most regions of France have relatively low deviations, with a few areas in the north and central parts showing higher deviations. The color scale on the right side of the map ranges from 0.6 to 1.8, indicating the magnitude of the deviations. The map highlights regions with significant excess mortality due to heat waves, particularly in the central and northern parts of the country.

To address this, we estimated separate regional mortality tables and computed the parameters $\alpha ^{(r)}_x$ , $\kappa ^{(r)}_t$ , $\beta ^{(r)}_x$ independently for each region $r \in \mathcal{R}$ . These regional parameters define the annual baseline hazard $\hat {\mu }^{(r)}_{x,t}$ , which serves as the foundation for the higher-resolution temporal model described in the following sections.

(6) \begin{equation} \hat {\mu }^{(r)}_{x,t} = \hat {m}^{(r)}_{x,t} = \exp \left ( \alpha _x^{(r)} + \beta _x^{(r)} \kappa _t^{(r)} \right ) \end{equation}

3.2 Specification of the seasonal component using a Fourier decomposition

Mortality data in France exhibit a pronounced seasonal pattern, with higher mortality rates observed during winter months. This phenomenon is primarily attributable to winter-related pathologies, such as influenza and pneumonia (Rau et al., Reference Rau, Bohk-Ewald, Muszynska and Vaupel2018).

Figure 3

The seasonality $s^{(climate\_type)}_{.,.,p}$ represented as a function of the time period $p$ .

Figure 3 Long description

The line graph displays seasonality as a function of time period for eight different climate types across various regions. The x-axis represents the 5-day period ranging from 0 to 72, while the y-axis indicates the seasonality values ranging from 0.85 to 1.20. Each line represents a different climate type, with distinct colors for each. The lines show fluctuations in seasonality over time, with notable peaks and troughs. All values are approximated.

Figure 4

The seasonality $s^{N}_{x,.,p}$ represented as a function of the time period $p$ .

Figure 4 Long description

The line graph displays seasonality as a function of the time period for various age groups. The x-axis represents the 5-day period ranging from 0 to 72, while the y-axis indicates seasonality values. The graph includes six distinct lines, each representing a different age group: 40 to 44, 50 to 54, 60 to 64, 70 to 74, 80 to 84, and 90 to 94. The lines show fluctuations in seasonality over time, with each age group exhibiting unique patterns. The seasonality values range from approximately 0.9 to 1.3. The lines intersect and diverge at various points, illustrating the varying impacts of seasonality across different age groups over the specified time periods. All values are approximated.

To calculate excess mortality at specific points throughout the year, seasonality must be formally accounted for. We define the expected number of deaths, $\widehat {d}^{\,(r)}_{x,t,p}$ , based on the exposure $E^{(r)}_{x,t,p}$ and the mortality rates $\hat {m}^{(r)}_{x, t}$ derived from the baseline model in Section 3.1:

(7) \begin{equation} \widehat {d}^{\,(r)}_{x,t,p} = E^{(r)}_{x,t,p} \left ( \frac {p \times \hat {m}^{(r)}_{x,t+1} + (72-p) \times \hat {m}^{(r)}_{x,t}}{72} \right ). \end{equation}

The empirical seasonality, $s^{(r)}_{x,.,p}$ is defined as the mean ratio of observed to expected deaths for each period p across the ten-year study period:

(8) \begin{equation} s^{(r)}_{x,.,p} = \frac {1}{\#\mathcal{T}} \sum _{t\in \mathcal{T}} \frac {d^{(r)}_{x,t,p}}{\widehat {d}^{\,(r)}_{x,t,p}} \end{equation}

The seasonal index $s$ thus depends on three parameters: age group $x$ , region $r$ , and time period $p$ . We examined these dependencies across regions and ages. First, seasonality was calculated for each of the eight climate types defined by Joly et al. (Reference Joly, Brossard, Cardot, Cavailhes, Hilal and Wavresky2010). As shown in Figure 3 (already mentioned in Section 2.2), seasonal patterns are remarkably consistent across all climate types. Second, seasonality was calculated by age group at the national level, denoted $s^N_{x,.,p}$ , as illustrated in Figure 4. The highest age groups have mortality rates that are more strongly dependent on seasonality. We can deduce from these two graphs that it is possible to construct a seasonality term based solely on age group and time step, and not on the region, calculated as:

(9) \begin{equation} s^{N}_{x,.,p} = \frac {1}{\#\mathcal{T}} \sum _{t\in \mathcal{T}} \frac {\sum _{r\in \mathcal{R}}d^{(r)}_{x,t,p}}{\sum _{r\in \mathcal{R}}\widehat {d}^{\,(r)}_{x,t,p}} \end{equation}

Following Robben et al. (Reference Robben, Antonio and Kleinow2025) and Schöley (Reference Schöley2021), we modeled the seasonality $s^{N}_{x,.,p}$ using a Fourier decomposition. The modeled seasonality is denoted $\tilde {s}^{N}_{x,.,p}$ . Our analysis indicates that while using only the first two harmonics is common, it is insufficiently complex for this dataset. Incorporating a third harmonic represents a statistically significant improvement in goodness of fit; however, the addition of a fourth or subsequent terms yields diminishing returns. Further details are provided in Appendix B. Consequently, we adopted a three-term Fourier series model, with coefficients calculated independently for each age group. By combining the annual regional trend with this seasonal component, we obtain the comprehensive baseline mortality term:

(10) \begin{equation} \hat {b}^{(r)}_{x,t,p} = \widehat {d}^{\,(r)}_{x,t,p} \times \tilde {s}^{N}_{x,.,p}. \end{equation}

This term $\hat {b}^{(r)}_{x,t,p}$ represents the expected death count under normal environmental conditions and is utilized as the offset in the integrated weather-mortality model (Section 3.3).

3.3 Modeling observed mortality with weather data

Justification for the incorporation of weather data: The mortality model is constructed in three sequential stages: the long-term trend (Section 3.1), the seasonal component (Section 3.2), and the weather-driven deviation from baseline.

To motivate the inclusion of meteorological covariates, we examine the standardized mortality ratio (SMR), which represents the relative deviation of observed deaths from the established baseline:

(11) \begin{equation} SMR_{x,t,p}^{(r)}\,:\!=\, \frac {d_{x,t,p}^{(r)}}{\hat {b}_{x,t,p}^{(r)}} \end{equation}

where $d_{x,t,p}^{(r)}$ denotes the observed number of deaths and $\hat {b}_{x,t,p}^{(r)}$ is the expected baseline mortality.

As illustrated in Figure 5, these ratios exhibit a clear dependency on temperature and humidity. The left-hand plot demonstrates that mortality remains relatively stable and close to the baseline ( $SMR\approx 1$ ) for temperatures between 12 $^\circ$ C and 19 $^\circ$ C. However, above 20 $^\circ$ C, the ratio begins to rise, transitioning into an exponential increase once temperatures exceed 26 $^\circ$ C. The right-hand plot indicates that, outside of the extremes (below 45% and above 93% relative humidity), the SMR is a decreasing function of humidity, characterized by a relatively stable, slightly negative slope.

Figure 5

Illustration of the relationship between the SMR and the weather data. Temperature and the left-hand side, and humidity on the right-hand side. This corresponds to the average of SMR relative to baseline $d^{(r)}_{x,t,p} / \hat {b}^{(r)}_{x,t,p}$ for every half of a degree celsius or for every percentage point of relative humidity.

Figure 5 Long description

A scatter plot with two panels. The left panel shows the relationship between the standardized mortality ratio and average mean temperatures over five-day intervals. The x-axis represents the average of the mean temperatures in degrees Celsius, ranging from negative ten to thirty. The y-axis represents the standardized mortality ratio relative to baseline, ranging from zero point seven to one point four. The right panel shows the relationship between the standardized mortality ratio and average mean humidity over five-day intervals. The x-axis represents the average of the mean humidity in percentage, ranging from thirty to one hundred. The y-axis represents the standardized mortality ratio relative to baseline, ranging from zero point seven to one point four. Each panel contains dozens of data points, showing varying trends and patterns. The left panel indicates a positive correlation between temperature and mortality ratio, while the right panel shows a less clear relationship between humidity and mortality ratio. All values are approximated.

Model specification: We assume that the observed death counts in region $r$ , for period $p$ and age group $x$ , are realizations of a Poisson-distributed random variable $D_{x,t,p}^{(r)}$ :

(12) \begin{align} D_{x,t,p}^{(r)} \sim \text{Poisson}\left ( \lambda _{x,t,p}^{(r)} \right ), \end{align}

where $\lambda _{x,t,p}^{(r)} = \mathbb{E}[D_{x,t,p}^{(r)}]$ denotes the expected number of deaths.

The objective is to account for the deviations between the baseline number of deaths, $\hat {b}^{(r)}_{x,t,p}$ , and the observed counts, $d_{x,t,p}^{(r)}$ , by incorporating the exogenous covariates defined in Table C.1. We assume a functional relationship where the baseline acts as an offset, such that:

(13) \begin{equation} \log \left ( \lambda _{x,t,p}^{(r)} \right ) = \log \left ( \hat {b}_{x,t,p}^{(r)} \right ) + h \left (W_{t,p}^{(r)}\right ) \end{equation}

where $W_{t,p}^{(r)}$ represents the vector of explanatory variables (including meteorological data) for region $r$ during period $p$ of year $t$ . Under this specification, the expected death counts are given by $\hat {\lambda }_{x,t,p}^{(r)} = \hat {b}_{x,t,p}^{(r)} \times \exp \left (\hat {h}(W_{t,p}^{(r)})\right )$ , where the multiplicative effect $\exp \left (\hat {h}(W_{t,p}^{(r)})\right ) - 1$ represents the percentage change in mortality attributable to external conditions. We estimate the unknown function $h({\cdot})$ using a machine learning framework.

The covariate vector $W_{t,p}^{(r)}$ encompasses a wide range of meteorological and pollution data. Our analysis is not restricted to temperature; consistent with previous literature (Kalkstein and Greene (Reference Kalkstein and Greene1997), Barreca (Reference Barreca2012)), we recognize that mortality is statistically dependent on broader climatic factors, such as humidity. However, to maintain model parsimony and predictive accuracy, we ultimately select a subset of variables focused on temperature and humidity. As discussed in Section 4.2, the inclusion of the exhaustive set of environmental indicators (listed in Appendix C) was found to degrade the performance of the machine learning algorithm due to increased noise and over-fitting.

4.1 Construction of the dataset for the machine learning process

Initial results, as illustrated in Figure 5, indicate complex, nonlinear relationships between meteorological variables and mortality rates. Furthermore, the inclusion of multiple weather indicators for each five-day period significantly increases the dimensionality of the covariate space. Consequently, we employ machine learning algorithms to capture these high-dimensional dependencies. Within our supervised learning framework, the target variable $y$ represents the observed death counts ( $y = d^{(r)}_{x,t,t}$ ), and the baseline mortality count ( $\hat {b}^{(r)}_{x,t,t}$ ) is used as an offset. In constructing the feature matrix $X$ , three primary considerations, detailed below, are addressed. The machine learning algorithm will estimate $\hat {y} = \hat {\lambda }^{(r)}_{x,t,p}$ .

Aggregation of daily data into five-d ay intervals . As the meteorological data are recorded daily, whereas the model operates on a five-day grid, temporal aggregation is required. For minimum and maximum temperatures, we retain the absolute maximum, minimum, and mean values over each five-day period. This approach identifies the full range of temperature fluctuations within each interval. Mean temperature and humidity indicators are aggregated using their arithmetic means. Given the relatively marginal impact of precipitation and wind on mortality rates reported in the literature, we include only two summary indicators for each of these variables. Additionally, we incorporate the wet bulb temperature (WBT), defined as the temperature of a parcel of air cooled to saturation (100% relative humidity) through the evaporation of water into it, with the latent heat supplied by the parcel itself (Dunlop (Reference Dunlop2008)). The full specification of these indicators is provided in Table C.1.

Lagged effects. To capture the delayed impact of weather on mortality during period $p$ of year $t$ , we include meteorological data from both the concurrent period and the immediately preceding period. This specification accounts for a window of up to ten days without excessively increasing the number of features. This choice is informed by research across 18 French cities, which suggests that the impact of heat is most acute within 0–3 days of exposure but can persist for up to 10 days (Pascal et al. (Reference Pascal, Wagner, Corso, Laaidi, Ung and Beaudeau2018)).

Harvesting effect. The “harvesting effect,” or mortality displacement, refers to a phenomenon where a temporary increase in death rates is followed by a period of lower-than-average mortality. This occurs when a severe event (such as a heat wave) primarily affects individuals who are already frail or nearing the end of their life. These individuals die earlier than they otherwise would have, causing a spike in deaths during period $p$ and a corresponding deficit in period $p+1$ . Accounting for this displacement is essential to ensure that the model does not overestimate the long-term impact of extreme weather events.

4.2 Machine learning for the observed deaths count modeling

As discussed previously, the relationship between weather variables and mortality is complex and inherently nonlinear. We, therefore, model the observed mortality using machine learning methods that can flexibly capture such relationships. In particular, the models incorporate both contemporaneous weather indicators and lagged indicators to account for delayed mortality effects. In addition, we explicitly include the SMR of the previous period $SMR^{(r)}_{x,t,p-1}$ as a predictor of $d^{(r)}_{x,t,p}$ , to capture the so-called “harvesting effect.”

Choice of models. We considered three types of models that can incorporate a Poisson specification to estimate $h$ with a function $\hat {h}$ of Equation 13:

• • a Poisson Generalized Linear Model (GLM),

• • XGBoost (Chen & Guestrin, Reference Chen and Guestrin2016), and

• • CatBoost (Prokhorenkova et al., Reference Prokhorenkova, Gusev, Vorobev, Dorogush and Gulin2018).

The GLM provides a transparent linear benchmark, while XGBoost and CatBoost are gradient boosting methods that allow for nonlinearities and interactions between predictors.

For the GLM, the statsmodel library (Seabold and Perktold (Reference Seabold and Perktold2010)) allows one to directly set the offset as a parameter of the GLM. The logarithms of the baseline death counts are used as an offset. Indeed, a Poisson GLM models the log of the expected count $\mu$ as a linear combination of predictors:

\begin{equation*} log(\mu )=X\beta + log(\textrm {offset})\end{equation*}

For both boosting algorithms, we adopted a log-link function, so that the logarithm of the expected mortality counts is modeled as the sum of the predictors and the log-offset derived from the baseline mortality. This ensures that the predictors act multiplicatively on the mean mortality count, which is essential for two reasons. First, scale invariance: the relative effect of predictors (e.g., a 10% increase in risk due to a heat wave) remains consistent regardless of whether the baseline mortality is low or high. Second, positivity: by construction, the exponential of the linear predictor guarantees that the expected counts are strictly positive, avoiding the theoretical issue of negative death counts that could arise in an additive formulation. In practice, this corresponds to specifying ‘objective’: ‘count:poisson’ and ‘eval_metric’: ‘poisson-nloglik’ for XGBoost, and ‘loss_function’: ‘Poisson’, ‘eval_metric’: ‘Poisson’ for CatBoost.

Hyperparameter tuning and validation. Model performance and generalization were evaluated using both in-sample (training) and out-of-sample (test) datasets. To determine the optimal hyperparameters for the machine learning framework, we employed a 5-fold cross-validation procedure on the training set using the Optuna optimization framework (Akiba et al. (Reference Akiba, Sano, Yanase, Ohta and Koyama2019)). Specifically, we utilized the tree-structured Parzen Estimator (TPE) (Watanabe (Reference Watanabe2025)), a Bayesian optimization approach, to minimize the average negative Poisson log-likelihood (PLL) across the validation folds. The search space included a broad range of parameters – such as the learning rate, tree depth, and various regularization terms, the details of which are provided in Appendix D. We adopted a “search before exploit” strategy by using the first thirty trials as exploration only. We also introduced a median pruner to terminate unpromising trials early, allowing for a deeper exploration of high-performing parameter combinations. Finally, the predictive accuracy of the models was compared against the benchmark baseline using three primary evaluation metrics: Mean Squared Error (MSE), PPLL, and Mean Poisson Deviance (MPD). The Normalized PLL (NPLL, PLL divided by the number of rows) was also computed to evaluate overfitting.

Variable s. Although gradient boosting methods are, in principle, capable of automatic feature selection, initial experiments with the full set of weather variables led to overfitting and degraded out-of-sample performance, even under intensive hyperparameter optimization with wide search ranges and large numbers of trials. This behavior is likely due to the strong collinearity among meteorological variables. We therefore adopted an iterative procedure in which features were added step by step and retained only if they improved performance consistently across folds. The final feature set comprised demographic variables (“age_group,” “sex”), temporal variables (“time_period”), geographical location (“latitude,” “longitude”), weather variables (“TM_mean,” “TM_mean_pm1,” “UM_mean,” “UM_mean_pm1”), and the harvesting effect proxy (“ecart_dc_tp”). All lagged features with the suffix “_pm1” refer to the values in the previous 5-day period. We restricted the main analysis to the summer period ( $p \in [30,54]$ ), when weather-related mortality effects are expected to be strongest. Weather data are detailed in C.

Table 3.

Comparative performance of the baseline (Offset) and machine learning models. Performance is evaluated using mean squared error (MSE), poisson log-likelihood (PLL), Normalized PLL (NPLL), and mean poisson deviance (MPD). The parameters used to obtain these results are detailed in D

Table 3. Long description

The table presents a comparative performance analysis of the baseline (Offset) and three machine learning models: GLM Poisson, XGBoost, and CatBoost. The performance metrics used are mean squared error (MSE), poisson log-likelihood (PLL), Normalized PLL (NPLL), and mean poisson deviance (MPD). The table is divided into two main sections: Training set and Test set, each containing five columns: MSE, PLL, NPLL, and MPD. There are four rows, each representing a different model. For the Offset model, the Training set shows MSE of 4.243, PLL of negative 563,119, NPLL of negative 1.778, and MPD of 1.087. The Test set for Offset shows MSE of 4.268, PLL of negative 140,716, NPLL of negative 1.777, and MPD of 1.088. For the GLM Poisson model, the Training set shows MSE of 4.174, PLL of negative 562,008, NPLL of negative 1.774, and MPD of 1.080. The Test set for GLM Poisson shows MSE of 4.196, PLL of negative 140,443, NPLL of negative 1.773, and MPD of 1.082. For the XGBoost model, the Training set shows MSE of 3.990, PLL of negative 559,310, NPLL of negative 1.765, and MPD of 1.063. The Test set for XGBoost shows MSE of 4.150, PLL of negative 140,302, NPLL of negative 1.771, and MPD of 1.078. For the CatBoost model, the Training set shows MSE of 4.137, PLL of negative 561,575, NPLL of negative 1.773, and MPD of 1.077. The Test set for CatBoost shows MSE of 4.174, PLL of negative 140,356, NPLL of negative 1.772, and MPD of 1.079. Notable trends include XGBoost showing the lowest MSE and NPLL in both training and test sets, indicating better performance compared to other models.

Figure 6

Comparative analysis of predicted and observed mortality. Predicted mortality (blue), observed mortality (orange), and the baseline mortality rate (green) are displayed alongside mean temperature (red, dotted). Data are aggregated across sex and all age groups, and for a specified region $(r=37)$ and year $(t=2015)$ .

Figure 6 Long description

The line graph presents a comparative analysis of predicted, observed, and baseline mortality rates over a specified time period. The x axis represents the time period, ranging from 30 to 55. The y axis on the left indicates the number of deaths, ranging from 600 to 1050. The y axis on the right shows the mean temperature in degrees Celsius, ranging from 12 to 28. The blue line represents the predicted number of deaths, the orange line shows the observed number of deaths, and the green line indicates the baseline number of deaths. Additionally, a red dotted line represents the mean temperature. The graph illustrates fluctuations in mortality rates and temperature trends over the given time period. All values are approximated.

4.3 Results and interpretation of the machine learning process

Table 3 summarizes the predictive performance of the baseline and machine learning models following hyperparameter optimization. The results indicate that XGBoost emerges as the most robust model, outperforming the alternative specifications across all evaluated metrics. Specifically, it achieves both the minimum MSE and the maximum PLL on both the training and test datasets. As evidenced by the results, the gradient boosting architectures (XGBoost and CatBoost) demonstrate superior predictive power compared to the GLM Poisson. While these complex models are often susceptible to overfitting, the narrow performance gap between our training and test sets suggests that the implemented regularization strategies were effective. To ensure strong generalization and mitigate the risk of over-fitting, we employed a rigorous Bayesian optimization framework for hyperparameter tuning, incorporating constraints on tree depth and $L1$ and $L2$ regularization alongside 5-fold cross-validation. Figure 6 demonstrates the sensitivity of the algorithm to temperature fluctuations. The significant spike in “TM_mean at period 36 directly correlates with the peaks in both observed and predicted deaths, highlighting the model’s capacity to capture extreme weather effects. The trajectories of temperature, observed mortality, and predicted deaths remain closely aligned throughout the remainder of the series, suggesting the model effectively captures the temperature fluctuations and their impact on mortality. Additional results regarding the model’s calibration relative to the baseline are presented in F for both the training and test sets.

Figure 7

Feature importance chart showing how much impact each feature has on the model output.

Figure 7 Long description

A scatter plot displays the impact of various features on a model’s output using SHAP values. The x-axis represents SHAP values, indicating the impact on the model output, while the y-axis lists the features. The features include sex, TM_mean, age_group, TM_mean_pm1, time_period, ecart_dc_tp, longitude, latitude, and UM_mean. Each feature is represented by a series of points, with colors ranging from blue to red indicating low to high feature values. The plot shows how each feature’s value influences the model’s predictions, with some features having a more significant impact than others. The data points are distributed along the x-axis, showing varying degrees of influence. The plot helps in understanding which features are most important in the model’s decision-making process. All values are approximated.

Figure 8

Shapley values given the values of the main weather conditions.

Figure 8 Long description

A scatter plot with two subplots labeled a and b. Subplot a shows the relationship between temperature mean in degrees Celsius and SHAP value for temperature mean. The x-axis represents temperature mean ranging from 5 to 30 degrees Celsius, and the y-axis represents SHAP value for temperature mean ranging from negative 0.05 to 0.15. The data points form a positive correlation, indicating that as the temperature mean increases, the SHAP value for temperature mean also increases. Subplot b shows the relationship between urbanization mean in percentage and SHAP value for urbanization mean. The x-axis represents urbanization mean ranging from 40 to 100 percent, and the y-axis represents SHAP value for urbanization mean ranging from negative 0.03 to 0.03. The data points form a less clear pattern with some clustering around certain urbanization mean values. All values are approximated.

Figure 9

Cross-effect of age (a) and latitude (b) with the mean temperature on the model, illustrated with the Shapley values.

Figure 9 Long description

A scatter plot illustrates the relationship between mean temperature in degrees Celsius and Shapley values, with age group and latitude as variables. The x-axis represents the mean temperature in degrees Celsius, ranging from 5 to 30 degrees. The y-axis represents Shapley values for the mean temperature, ranging from negative 0.05 to 0.15. The data points are color-coded, with red representing higher age groups and blue representing lower age groups. The left plot shows the cross-effect of age group on the model, while the right plot shows the cross-effect of latitude. Both plots exhibit a positive correlation between mean temperature and Shapley values. The data points form clusters and show a general upward trend as the mean temperature increases. All values are approximated.

The marginal contributions of each feature to the model’s output were quantified using Shapley values, implemented via the SHAP library (Lundberg and Lee (Reference Lundberg and Lee2017)). Derived from cooperative game theory, these values provide a mathematically rigorous framework for attributing the contribution of individual “players” (features) to the final result. Analysis of the feature importance hierarchy (Figure 7) reveals that gender is the primary discriminating variable, exhibiting a distinct separation between male (blue) and female (red) cohorts. This observation aligns with existing literature regarding sex-based physiological sensitivities to thermal stress (Folkerts et al. (Reference Folkerts, Bröde, Botzen, Martinius, Gerrett, Harmsen and Daanen2022)). The influence of temperature on mortality is characterized by two distinct phases. Below a 15 $^\circ$ C threshold, the SHAP values remain stable, indicating a negligible impact on excess mortality. Conversely, above 15 $^\circ$ C, the values become a strictly increasing function of temperature, demonstrating an accelerating, almost exponential, effect on mortality risk (Figure 8).

While the SHAP summary plot provides a global overview, it can be laterally reductive. For instance, Figure 8 might erroneously suggest that younger individuals are more vulnerable, as the “age group” summary depicts red dots to the left of blue dots. This discrepancy arises because summary plots aggregate average marginal contributions and often obscure complex interaction effects. To resolve this, dependence plots were employed to isolate cross-effects (Figure 9). These plots confirm that age acts as a significant force multiplier; at elevated temperatures, the risk is concentrated among older individuals (represented by the vertical divergence of red dots in Figure 9 a). A similar interaction is observed with latitude (Figure 9 b), where higher latitudes correlate with increased excess mortality during heat events. This suggests a lack of thermal resilience compared to populations in Southern France, who benefit from more frequent acclimatization.

Finally, the interaction between temperature and humidity (Figure 10) underscores the detrimental nature of humid heat. While France’s temperate climate typically precludes the simultaneous occurrence of extreme heat and high humidity, the Mediterranean sub-analysis (Figure 10 b) reveals a cluster of high-humidity, high-temperature points corresponding to the highest observed Shapley values. These findings necessitate targeted public health interventions in Mediterranean regions, where the synergistic effect of humidity and heat poses the most acute mortality risk.

Figure 10

Cross-effects of temperature and humidity in France (a) compared with the area around the Mediterranean see (b).

Figure 10 Long description

A scatter plot showing the relationship between temperature and humidity in France and the Mediterranean region. The plot consists of hundreds of data points, with the x-axis representing the mean temperature in degrees Celsius and the y-axis representing the mean humidity percentage. The data points are color-coded, with blue and pink dots indicating different datasets or conditions. The left plot focuses on France, with temperature values ranging from 5 to 30 degrees Celsius and humidity values ranging from 55 to 85 percentage. The right plot focuses on the Mediterranean region, with temperature values ranging from 12.5 to 30 degrees Celsius and humidity values ranging from 55 to 75 percentage. Both plots show a positive correlation between temperature and humidity, with data points clustering more densely as temperature increases. The scatter plots highlight the cross-effects of temperature and humidity in these regions, providing insights into how these environmental factors interact. All values are approximated.

5.1 Method 1: Worst scenario using lots of climate projections

The first methodology employs an ensemble of regional climate projections. Specifically, twelve scenarios corresponding to the RCP 8.5 pathway were obtained from the DRIAS website.Footnote ² Developed by Météo-France in collaboration with the French climate science community (including IPSL, CERFACS, and CNRM), DRIAS provides high-resolution, regionalized climate data tailored for adaptation strategies within France.

By combining 15 years of historical data with these twelve projections, a synthetic pool of 180 distinct summer seasons was generated. From this ensemble, the summer exhibiting the highest mean temperature was isolated for further analysis. This extreme case is derived from a projection utilizing the HadGEM2 global circulation model (Collins et al. (Reference Collins, Bellouin, Doutriaux-Boucher, Gedney, Hinton, Jones, Liddicoat, Martin, O’Connor and Rae2008)), with downscaling performed via the CCLM4-8-17 regional model.

The primary advantage of utilizing these raw climate scenarios lies in the inherent physical consistency between variables. Because these projections are underpinned by fundamental geophysical equations and rigorous statistical regionalization, they ensure that the simulated extreme events, while statistically unlikely, remain meteorologically plausible. Consequently, the direct application of these scenarios provides a robust guarantee of physical coherence across the model features.

5.2 Method 2: Worst scenario using the most extreme possible temperatures

To simulate summers with extreme temperatures, the extreme values theory can be useful (Coles et al., Reference Coles, Bawa, Trenner and Dorazio2001). In fact, historical data are not numerous enough and climate models have difficulties in representing the tails of the distributions (Van Oldenborgh et al. (Reference Van Oldenborgh, Wehner, Vautard, Otto, Seneviratne, Stott, Hegerl, Philip and Kew2022), Vautard et al. (Reference Vautard, Cattiaux, Happé, Singh, Bonnet, Cassou, Coumou, D’Andrea, Faranda, Fischer, Ribes, Sippel and Yiou2023)). Extreme value theory has been used to simulate extreme temperatures (Wehner et al. (Reference Wehner, Stone, Shiogama, Wolski, Ciavarella, Christidis and Krishnan2018), Auld et al. (Reference Auld, Hegerl and Papastathopoulos2023), Van Oldenborgh et al. (Reference Van Oldenborgh, Wehner, Vautard, Otto, Seneviratne, Stott, Hegerl, Philip and Kew2022)). However, simply applying it in a classical framework is not adequate. In the context of temperatures, the extreme value theory has indeed many shortcomings:

• • Using the Block-Maxima approach (BM) requires many years of historical data, mixing data from different climates, but the distribution is not stationary and changes with global warming. The smoothed global mean surface temperature (GMST) is often used to account for non-stationarity (Wehner et al. (Reference Wehner, Stone, Shiogama, Wolski, Ciavarella, Christidis and Krishnan2018), Philip et al. (Reference Philip, Kew, van Oldenborgh, Otto, Vautard, van der Wiel, King, Lott, Arrighi, Singh and van Aalst2020)). The location parameter of the generalized extreme value distribution (GEV) is then written:

  (14) \begin{equation} \mu = \mu _0 + \mu _1 * GMST \end{equation}
  The scale and shape parameters are often assumed constant.

• • Using the Peak-Over-Threshold method (PoT) gives uncertain results. In general, there are no clear thresholds that emerge from the analysis of the mean residual life and the stability of the parameters, as illustrated in Figure 11. In the mean residual plot 11 c), the uncertainty around the estimation is important since the first threshold. In the stability plot of parameters 11 a) the shape parameter $\xi$ oscillates very slightly around 0, which greatly changes the obtained distribution. And the shape parameter $\sigma$ , which must be negative, tends toward 0 as it can be seen in Figure 11 b).

• • Once the data to fit the GEV have been chosen by the PoT or BM approach, the shape parameter of the GEV distribution is almost always found to be negative in heat-wave analyses. This means that the resulting distribution has an upper bound. But, as written in Noyelle et al. (Reference Noyelle, Robin, Naveau, Yiou and Faranda2024): “Recent intense heat waves […] have nonetheless challenged the reliability of the estimation of the statistical upper bound by breaking it sometimes by a large margin.” To alleviate this shortcoming, Zhang and Boos (Reference Zhang and Boos2023), Noyelle et al. (Reference Noyelle, Zhang, Yiou and Faranda2023), and Noyelle et al. (Reference Noyelle, Robin, Naveau, Yiou and Faranda2024) have studied the physics of heat waves and incorporated it into the GEV model. We use their work to generate data that is specifically aligned with our research objectives.

Figure 11

Study of the threshold for the POT method on historical data measured at the Paris Montsouris station with the pyextreme library. It shows both the estimation (red line) and the confidence interval at 95%.

Figure 11 Long description

The line graph presents the study of the threshold for the POT method on historical data measured at the Paris Montsouris station using the pyextreme library. The x axis represents the threshold values ranging from 26 to 38. The y axis in the first graph represents the shape values ranging from negative 2 to 4, and in the second graph, it represents the modified scale values ranging from negative 100 to 100. The third graph shows the mean excess values ranging from 0.5 to 3. The red line indicates the estimation, and the blue shaded area represents the confidence interval at 95%. All values are approximated.

In Zhang and Boos (Reference Zhang and Boos2023), the authors highlighted that high temperatures are physically limited by the moist convective instability of the air column. The convective instability triggers precipitation when the temperature reaches this limitation, and the temperature falls down due to precipitation. This physical limit defines a relevant upper bound for surface temperatures. This was observed in July 2019 in the Paris region, for example. The temperature gradually rose from July 21 to July 25 until it reached more than 40 $^\circ$ C on July 25. The following day, it began to rain and the maximum temperature dropped by 10 $^\circ$ C. Noyelle et al. (Reference Noyelle, Robin, Naveau, Yiou and Faranda2024) propose to estimate statistical GEV models on time series of yearly maxima where the upper bound of the distribution is imposed by the physical knowledge of the system. The physical upper bound can be calculated using the equation:

(15) \begin{equation} T_{s, max} = T_{500} + \frac {L_v}{c_p} ( Q_{sat}(T_{500}) - Q_s) + \frac {g}{c_p}(Z_{500} - Z_s) \end{equation}

where $T_{s, max}$ is the maximal reachable temperature at the surface( $s$ ), $T_{500}$ is the air temperature at 500 hPa, $L_v$ is the latent heat of vaporization, $c_p$ is the specific heat of air at constant pressure, $Q_{sat}(T_{500})$ the saturation specific humidity at $T_{500}$ , $Q_s$ the surface specific humidity, $g$ the gravitational constant, $Z_{500}$ the geopotential height at 500 hPa and finally $Z_s$ the elevation of the surface. The data used to compute $T_{s, max}$ come from Copernicus European Regional ReAnalysis (CERRA) datasets, more specifically the sub-daily reanalysis data for Europe on pressure levels from 1984 to present.Footnote ³ As there was no open source available data for projections that could be easily be founded including these variables, we use the reanalysis data from 2010 to 2019 to compute the upper bounds. The physical upper bound was calculated for each day, and we retain the maximum of these bounds as $T_{s, max}$ . The results are shown in Figure 12. 50 $^\circ$ C may be reached at various locations in France, mainly in the south of the country. The Paris region could see peaks of 46 $^\circ$ C. These are maximums calculated in the current climate. With global warming, these upper limits could rise further.

Figure 12

$T_{s, max}$ for each cluster computed on the reanalysis data.

Figure 12 Long description

A heat map of France showing temperature distribution with a color scale ranging from 35 to 55 degrees Celsius. The map uses a gradient color scheme where lighter colors represent lower temperatures and darker colors represent higher temperatures. The highest temperatures, indicated by dark red and maroon colors, are concentrated in the southern and central regions of France. The northern and western regions exhibit cooler temperatures, shown in lighter yellow and orange shades. The map provides a visual representation of temperature variations across different regions of France, highlighting areas with significantly higher or lower temperatures.

$T_{s, max}$ correspond to the upper bound for $TX\_max$ (cf. Appendix C for the notations) within the framework of our study. But the model finally uses only $TM\_mean$ , $TM\_mean\_pm1$ , and $UM\_mean$ as weather data. We calculated the difference between $TX\_max$ and $TM\_mean$ for each of the clusters for the hottest 5% of days during the summer period. This difference allowed us to recalculate the $TM\_mean$ values for the $T_{s, max}$ values. The method involves replacing the hottest 5-day period of the summer with newly calculated data from $T_{s, max}$ . This applies to all projection data.

5.3 Method 3: Stochastic weather generator for extreme heat wave

The third method consists of projections with longer-lasting heat. In Yiou et al. (Reference Yiou, Cadiou, Faranda, Jézéquel, Malhomme, Miloshevich, Noyelle, Pons, Robin and Vrac2023), the authors try to anticipate the worst heat wave that could have occurred during the Paris 2024 Olympics. They use the method described in Yiou and Jézéquel (Reference Yiou and Jézéquel2020). An importance sampling method is employed to adapt a stochastic weather generator for the projection of particularly long-lasting heat waves. A stochastic weather generator is a statistical model that produces realistic physical characteristics while requiring minimal computational resources (Ailliot et al., Reference Ailliot, Allard, Monbet and Naveau2015). Heat waves are characterized by persistent episodes of high temperatures that are associated with persistent anticyclonic atmospheric patterns. The objective is to construct coherent atmospheric circulation scenarios that give rise to these persistent anticyclonic characteristics. This corresponds to a Markov chain of temperature, with latent states of the atmospheric circulation. We used the temperature $T$ and the geopotential height at 500 hPa ( $Z_{500}$ ). A data set $\mathcal{E} = (\mathcal{Z}, \mathcal{T})$ was created from several models over several years where:

• • $\mathcal{Z}$ is the set of $Z_{500}^{t,m}$ data. $Z_{500}^{t,m} = \{Z_{500,\text{lat},\text{lon}}^{t,m}\}$ , we consider the surface of all points for all latitudes $lat$ and longitudes $lon$ of the geographical perimeter illustrated in Figure 13.

• • $\mathcal{T}$ the set of temperatures $T^{t,m}$ for time $t$ for the data issued by model $m$ . $T^{t,m}$ are the spatial averages of the mean daily temperature over hexagonal France.

Figure 13

Example of the field of the $Z_{500}^{t,m} = \{Z_{500,\text{lat},\text{lon}}^{t,m}\}$ for July 1, 2026, in the scenario from the IPSL-IPSL-CM5A-MR model.

Figure 13 Long description

A heat map titled ’Map of Z_500 values’ displays the distribution of Z_500 values across Europe and surrounding regions. The map uses a color gradient ranging from purple to yellow, indicating values from 5400 to 5900. The x-axis and y-axis are not explicitly labeled but represent geographical coordinates. Higher values are shown in yellow, while lower values are depicted in purple. The map highlights a gradient where the northern regions, particularly Scandinavia, exhibit lower values in purple, transitioning to higher values in green and yellow in central and southern Europe. The heat map provides a visual representation of atmospheric conditions, with notable variations across different regions.

The procedure adapted to our specific case is described in pseudocode in the algorithm in Figure 14. The projection for July is given as an example. One projection was made per month from June to September. The parameter that controls the weights toward high temperatures is $\alpha$ . An increase in $\alpha$ value is indicative of a greater significance being attributed to elevated temperatures.

Figure 14

Algorithm used for the synthetic weather generator.

Figure 14 Long description

The equation presents a synthetic weather generator with specific inputs and initialization steps. The decay rate alpha is set to 0.5, and the number of analogues Z 500 fields K is 20. The initialization involves calculating S K using an exponential function with alpha and K, followed by determining the normalization constant A. Weights are then calculated using A and an exponential function with alpha and a set of integers from 1 to K. The number of days in July, T, is set to 31, and a starting day in July is randomly selected. The main loop involves selecting the K closest analogues based on Euclidean distance, ranking them by temperature, sampling the next day using the weights, and updating the current day.

This method makes it possible to generate permanently high-temperature projections while maintaining consistency in atmospheric circulation.

5.4 Number of deaths under each of the extreme projection

The future meteorological trajectories generated by the three projection methodologies (illustrated in Figures 15 and 16) were employed as inputs for the model described in Section 4. These were used to quantify the excess mortality associated with each projection for the study population.

Figure 15

Comparison of the average TM_mean over the whole period ( $p\in [30;54]$ ), between the historical year for which this average is highest (2018) and the averages obtained by the projection methods (5.1 and 5.3).

Figure 15 Long description

A heat map compares the average TM_mean values across France for the historical year 2018 and two projection methods. The map is divided into three sections: Historical data (2018), Method 1, and Method 3. Each section shows a color gradient representing the average TM_mean values, with a scale ranging from 10 to 40. The historical data section shows predominantly lighter colors, indicating lower average values. Method 1 shows a mix of light to medium colors, suggesting a moderate increase in values. Method 3 displays predominantly darker colors, indicating higher average values. The color intensity increases from left to right, reflecting higher TM_mean values in the projection methods compared to the historical data.

Figure 16

Comparison of the maximum TM_mean over the whole period ( $p\in [30;54]$ ), between the historical year for which this maximum is highest (2019) and the maximums obtained by the projection methods (5.1, 5.2, and 5.3).

Figure 16 Long description

The heat map consists of four panels, each representing a different dataset. The top left panel shows historical data from 2019, with a color scale ranging from 10 to 40. The top right, bottom left, and bottom right panels display projections from Method 1, Method 2, and Method 3, respectively. Each panel uses a similar color scale to indicate the intensity of the TM_mean values across different regions of France. The color gradient moves from lighter shades, indicating lower values, to darker shades, indicating higher values. The historical data shows a varied distribution of TM_mean values, with some regions exhibiting higher intensities. The projection methods show different patterns of intensity distribution, with Method 1 and Method 2 displaying more concentrated areas of high values, while Method 3 shows a more dispersed pattern. The heat map provides a visual comparison of how the TM_mean values are projected to change over the period, highlighting regions with significant variations.

As demonstrated in Figure 15, mean temperatures under Method 1 exceed historical observations, with Method 3 yielding even higher elevations. Specifically, the national average for France was 19.5 $^\circ$ C in 2018, compared to 22.5 $^\circ$ C under Method 1 and 25.5 $^\circ$ C under Method 3. These variances, while seemingly modest, are highly significant given the nonlinear relationship between temperature and mortality previously identified in Figure 8. Results for Method 2 are omitted from Figure 15 as they are restricted to a specific five-day window.

Conversely, Figure 16 shows that the highest peak temperatures are, as anticipated, attained via Method 2, which represents the maximum theoretical thresholds. Methods 1 and 3 approach these limits, with all three methodologies simulating unprecedented temperature maxima. This underscores the methodological advantage of incorporating synthetic climate projections rather than relying solely on historical datasets. The observed congruence between Methods 1 and 3 is consistent with expectations; although derived from different datasets, both are produced by the same class of climate model, with Method 3 representing a specific stochastic draw from the climate model ensemble.

Table 4.

Mortality shocks in each scenario based on the three methods to generate hot summers

Table 4. Long description

The table presents data on mortality shocks under different scenarios based on three methods to generate hot summers. It includes two key metrics: heat-related excess mortality on the scope of the study and annualized heat-related excess mortality on the whole population. The table has four columns: Stress test ACPR, Method 1, Method 3, and Method 2 + 3. Row 1: Heat-related excess mortality (on the scope of the study), Stress test ACPR, 4.01 percentage, Method 1, 8.67 percentage, Method 3, 8.69 percentage, Method 2 + 3. Row 2: Annualized heat-related excess mortality on the whole population, Stress test ACPR, 0.43 percentage, Method 1, 1.02 percentage, Method 3, 2.20 percentage, Method 2 + 3, 2.21 percentage.

The sensitivity tests were executed using the newly generated synthetic meteorological datasets, maintaining the geographic and demographic distributions observed in 2019. Four distinct simulations were conducted.

In the primary validation test, meteorological data for each period $p$ were derived by averaging observations across the historical horizon (2010–2019). This resulted in a negligible deviation of 0.1% from the baseline mortality, providing further empirical validation of the model’s stability. Subsequently, excess mortality was projected for three scenarios: Method 1, Method 3, and a hybrid approach combining Methods 2 and 3. In this latter configuration, the highest-temperature week within the Method 3 series was substituted with the extreme values derived from Method 2.

Excess mortality relative to the baseline was calculated across all defined age cohorts and time periods. To facilitate a comparison with the shock proposed in the ACPR 2023 climate stress test, these results were annualized against the general population.

The findings, summarized in Table 4, indicate that the projected shocks are significantly more severe than those currently proposed by regulatory authorities. Indeed, the modeled shocks are up to five times greater than the ACPR’s projected levels. Analysis suggests that the sustained periods of elevated temperature projected in Method 3 exert the most profound influence on mortality rates. Conversely, the integration of Method 2 yielded a marginal impact (0.01%), likely due to the limited five-day window in which the temperatures were adjusted upwards.

The ACPR’s projections are predicated on the observed consequences of the 2022 heat waves, under the implicit assumption that future climatic events will not exceed the severity of those already recorded. In the context of accelerating anthropogenic climate change, such an assumption remains questionable. Furthermore, the shock proposed by the ACPR appears notably conservative when compared to findings from Santé Publique France (Pascal et al. (Reference Pascal, Wagner and Lagarrigue2023)). Their research indicates that between June 1 and September 15, 2022, 4.1% of all-cause mortality was attributable to heat, a 1.91% increase relative to the 2014–2021 average. When annualized, this increase corresponds to a mortality shock of approximately 0.48%, which is 12% higher than the benchmark proposed by the ACPR. By simulating unprecedented meteorological conditions through the diverse projection methodologies employed in this study, we are able to construct more prudent scenarios. Such an approach is arguably better suited to quantifying the tail-risk probabilities and the upper extreme of the mortality risk distribution, providing a more robust framework for solvency requirements.