Quantifying cumulative heat effects on heat-related mortality using MODIS land surface temperature data

2.1 Mortality data and study area

We obtained weekly all-cause mortality data stratified by gender and age group for the warm season (weeks 18–40) from the European Union Statistical Office (https://ec.europa.eu/) for the period from 2000 to 2023. The cumulative heat effects may not be immediately observable and can persist for weeks. Therefore, we used weekly aggregated mortality data to reduce short-term noise (e.g. weekend effects and random daily fluctuations) and better capture these multiweek impacts. The final dataset covers 26 European countries and 122 NUTS-3-level administrative units (basic information can be found in Supplementary Material Table S1), with a total of 7,289,587 death records (3,636,604 women, 3,652,983 men, no personal privacy involved), representing a population of over 164 million (Figure 1a). The study area includes NUTS-3 level administrative units from Albania (2), Austria (3), Belgium (6), Bulgaria (3), Czech Republic (3), Estonia (2), Finland (3), France (12), Greece (3), Hungary (5), Italy (10), Latvia (4), Lithuania (3), Netherlands (8), Norway (4), Poland (9), Portugal (5), Romania (5), Slovakia (4), Spain (12), Sweden (3), Switzerland (3), and the United Kingdom (7), with the sample size for each region indicated in parentheses. These regions are all densely populated and economically developed major cities/regions. We also selected an additional 90 cities as independent verifications to validate the model's generalization ability, representing a population of over 30 million. The specific spatial distribution is shown in Supplementary Material Fig. S1. This approach of using weekly mortality counts is supported by recent methodological research, which demonstrates its utility in obtaining stable estimates of temperature-mortality associations (Ballester et al. Citation2023; Wu et al. Citation2025).

Figure 1. Study region and data. (a) Spatial distribution of 122 major regions/cities in Europe, with GIS data for administrative boundaries obtained from Eurostat. (b) National climates and seasons, based on satellite-derived LST data from 2000 to 2023.

Read detailed figure description

The two panel diagram shows European regions with mortality data and climate season with LST data from 2000 to 2023. Panel a displays a map of Europe with major regions and cities. A color scale indicates mortality in deaths per million, ranging from less than 100 (light blue) to greater than 250 (dark red). An inset map shows Europe within a global context. A horizontal scale bar indicates kilometers from 0 to 1200. Panel b shows a vertical bar graph of climate season in Europe. Months from January to November are listed on the vertical axis. A red bar highlights the warm season from May to September, corresponding to weeks 18 to 40. To the right, two stacks of maps represent 504 daytime and 504 nighttime images of Europe, illustrating 2000 to 2023 LST data from MODIS. Below these, two maps of Europe show pixel level maximum LST. The top map shows daytime LST ranging from 12.7 degrees Celsius to 58.3 degrees Celsius. The bottom map shows nighttime LST ranging from 2.3 degrees Celsius to 34.5 degrees Celsius.

Display full size

During data processing, there were no missing mortality records within the available time range (2000–2023) for the cities included in this study. However, some cities had gaps in coverage for certain years within this period. No data imputation was applied, as the mortality data were complete for the time spans covered. To ensure clarity, we have listed the range of mortality data covering all 212 cities in this study in Supplementary Material Table S2. The mortality data used in this study were obtained from official public health or statistical agencies. All data were fully deidentified and anonymized before being accessed by the researchers, containing no personally identifiable information. The analysis was conducted at the aggregated population level, preventing any individual-level identification. As the study relies solely on the secondary use of anonymized public data without direct human participation, no ethical approval or informed consent was required. Data handling complied with the privacy and confidentiality regulations of the corresponding data providers and was used exclusively for scientific research purposes.

2.2 Temperature data

This study utilizes the MOD11A1 dataset (https://ladsweb.modaps.eosdis.nasa.gov/) from NASA's Terra satellite (1 km spatial resolution, 2000–2023 time series). Daily observations of daytime land surface temperature (D_LST) and nighttime land surface temperature (N_LST) were retrieved through thermal infrared channels. Due to cloud cover, satellite-derived LST observations may contain missing values, potentially introducing uncertainty in temperature exposure assessment. To evaluate the robustness of MODIS LST-based estimates and to assess whether the main findings were sensitive to cloud-related data gaps, ERA5-Land was incorporated as an independent, spatially continuous reference dataset. The two datasets were used in parallel within the same modeling framework, with the comparison focused on consistency in model performance and inferred heat–mortality relationships rather than direct temperature equivalence or pixel-by-pixel agreement. This parallel comparison provides evidence that the main conclusions are not sensitive to cloud-related data gaps in MODIS LST products. This design also allows us to examine whether satellite-derived LST can provide added value for large-scale heat-health assessment beyond a widely used reanalysis product. We also provided the results of linear fitting between MODIS and ERA5 at the site level (Supplementary Material Tables S3 and S4). A detailed description of the handling of cloud-affected pixels and the parallel validation strategy has been added to the Materials and Methods section.

ERA5-Land temperature data were obtained from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/) and provide hourly 2 m near-surface air temperature fields at an approximate spatial resolution of 9 km. Differences in spatial resolution between MODIS LST and ERA5-Land were addressed by extracting all available pixels or grid cells within each city or regional boundary and constructing exposure metrics at the aggregated spatial-unit level. For the MOD11A1 product, which provides one daytime and one nighttime observation per day, weekly maximum daytime (D_LST) and nighttime (N_LST) temperatures were derived by selecting the highest observed values within each calendar week (Figure 1b). For ERA5-Land, hourly temperatures were first converted to daily maximum daytime and nighttime values for each grid cell and then aggregated over a 7-day moving window to derive weekly maximum temperatures, ensuring consistency with the MODIS-based exposure metrics. Regional or city-level weekly maximum temperatures were obtained by averaging across all pixels or grid cells within the corresponding spatial units. Because ERA5-Land and MODIS characterize different thermal quantities—near-surface air temperature versus land surface temperature—the ERA5-based analysis was used as a reference comparison to evaluate the robustness of cumulative heat–mortality associations across different temperature representations, rather than as a direct substitute for LST.

The use of weekly temperature and mortality data reflects a deliberate methodological choice aimed at capturing the cumulative and delayed health impacts of sustained heat exposure, while reducing short-term noise caused by day-to-day variability, reporting irregularities, and missing observations. Recent methodological studies have demonstrated that, under appropriate modeling frameworks, weekly aggregation produces unbiased and consistent estimates of temperature-related relative risks, with performance comparable to daily models, differing mainly in slightly wider confidence intervals (Ballester et al. Citation2024; Basagaña and Ballester Citation2024). This temporal scale is therefore well suited for analyzes focusing on prolonged thermal stress rather than solely on acute heat extremes. During data preprocessing, the NASA LAADS DAAC platform was used to ensure data quality; cloud-contaminated pixels were removed using the Fmask algorithm (Zhu, Wang, and Woodcock Citation2015), outliers were filtered using the standard deviation method, and residual noise was reduced using the Savitzky–Golay smoothing filter (Huang et al. Citation2021). The death and temperature data used in this study are available at Qian (Citation2025).

2.3 Accumulated high temperature and modeling

To investigate the impact of high temperatures on mortality rates in major European cities, we employed a quasi-Poisson regression model incorporating nonlinear functions to estimate the exposure-response relationship between temperature and mortality across cities. The minimum mortality temperature (MMT), defined as the temperature corresponding to the lowest relative risk of mortality, was identified as the point of minimal risk within this function. Subsequently, the relative risk (RR) at each temperature was calculated as the ratio of the predicted mortality risk at that temperature to the risk at the MMT (Gasparrini et al. Citation2015; Gasparrini, Armstrong, and Kenward Citation2010). The detailed computational workflow for RR estimation is provided in the Supplementary Materials. To assess the cumulative heat effects on mortality, this study employed a machine learning-based modeling approach (detailed in Section 2.4) that simultaneously accounts for both the current week's temperature and temperature deviations over the preceding three weeks. Temperature deviations were calculated as the difference between each city's weekly average temperature and its MMT. And the DLNM-based RR analysis was used primarily for epidemiological interpretation of exposure–response patterns, whereas the RF model was designed to quantify multiweek cumulative heat contributions and to evaluate predictive performance across heterogeneous regions.

The current week's daytime maximum temperature and nighttime maximum temperature are represented as D_LST and N_LST. In contrast, the previous week's temperatures are D_LST_1 and N_LST_1, so the daytime cumulative high temperature of last week is D_LST_AT1 (D_LST_1 minus MMT_D), the nighttime cumulative high temperature of last week is N_LST_AT1 (N_LST_1 minus MMT_N). Similarly, the temperatures from two and three weeks ago are represented as D_LST_AT2, N_LST_AT2, and D_LST_AT3, N_LST_AT3, respectively. We defined the cumulative temperature deviation over four weeks as a combination of these temperature variables, representing the cumulative contribution of multiweek temperature exposure to the health risk for the current week.

In the machine learning model, we did not directly input the past weeks' temperature deviations but instead introduced a decay factor to adjust the intensity of their impact on mortality risk. Here, X*denotes the decay coefficient (0 ≤ X ≤ 1) applied to temperature deviations from previous weeks, representing the proportion of accumulated temperature (PAT) retained from earlier exposures. The key assumption is that more recent temperature exposures exert a stronger influence on mortality than older exposures. Accordingly, the decay factor decreases progressively over time, such that temperature deviations from the current week receive full weight (100%), while those from earlier weeks are downweighted toward zero (input variables are listed in Table 1). Furthermore, correlation analysis indicates that there are no significant multicollinearity issues among the input variables; therefore, all variables were included in the machine learning model (Supplementary Material Table S5). The process is as follows (using daytime as an example, and the current weekly average temperature and temperature difference are always used as inputs):

2.4 Simulation analysis and scenario design for cumulative heat effects

To further evaluate the cumulative heat effects and the role of heat-event duration, we conducted simulation analyzes based on the trained RF model. First, a counterfactual substitution experiment was performed by replacing the weekly high-temperature series (D_LST and N_LST) in 2022 with the corresponding weeks from 2019, and then predicting mortality using the RF model. Because accumulated temperature predictors require antecedent information, the temperature data of the first three weeks of 2019 were retained as initial lag inputs to ensure consistent calculation of accumulated temperature variables. The predicted mortality distribution after substitution was compared with the observed mortality in 2019 and 2022 to quantify the contribution of extreme heat conditions.

In addition, four synthetic exposure scenarios (A–D) were constructed to examine how the temporal continuity of heat events influences cumulative effects. These scenarios were designed with comparable heat intensity but different temporal patterns, including continuous heat exposure, intermittent exposure, and alternating/sporadic heat events with recovery intervals. For each scenario, the corresponding temperature series were used as model inputs, accumulated temperature variables were recalculated accordingly, and the trained RF model was applied to estimate weekly mortality. The average predicted mortality during heat-affected weeks was used to compare the cumulative impacts across scenarios.

2.5 Machine learning

Random Forest (RF) is an ensemble learning framework developed by Breiman, which improves prediction performance by constructing heterogeneous decision tree groups through bootstrap aggregation (Breiman Citation2001; Qian et al. Citation2024a). The algorithm generates multiple parallel training subsets through resampling with replacement and randomly selects a subset of features at each node split to enhance model diversity. This ensemble structure addresses the overfitting problem of single decision trees and integrates the predictions through majority voting (for classification tasks) or mean regression (for regression tasks), forming a robust model with strong generalization capabilities. In this study, hyperparameters such as the forest size, tree depth, and minimum sample size for leaf nodes were optimized to ensure the model's performance and computational efficiency. These hyperparameters were selected using an automated tuning approach to balance model complexity and predictive accuracy. The performance of the tuned model was evaluated based on metrics like accuracy and out-of-bag error to ensure robustness. Permutation importance method is employed for feature evaluation, which quantifies each feature's contribution to the model's predictive power by assessing the change in prediction performance when each feature is perturbed (Aldrich Citation2020; Qian et al. Citation2024b). The details of other comparison models (Poisson regression and XGBoost), hyperparameter tuning, and permutation importance can be seen in the Supplementary Materials.

Recent studies have demonstrated the feasibility of using RF to model weekly mortality data and to capture nonlinear and cumulative temperature–mortality relationships (Boudreault et al. Citation2024; Kim et al. Citation2019; Robben, Antonio, and Kleinow Citation2025). To isolate short-term temperature variability from multiweek trends and nonclimatic disturbances, we implemented a structured comparison framework rather than relying on a single model specification. Specifically, time-only models were used to characterize temporal patterns unrelated to temperature, temperature-only models to quantify temperature-driven variability, and combined models to evaluate the incremental explanatory contribution of temperature beyond temporal structure. Temporal shuffling tests were further conducted as sensitivity analyzes to confirm that observed associations were not artifacts of temporal autocorrelation.

To account for potential nonclimatic structural changes, including long-term warming trends and pandemic-related mortality shocks, an additional pre-COVID model (2000–2019) was trained and compared with the full-period model. Age-stratified models were developed as a demographic sensitivity analysis to account for population aging and age-specific vulnerability, rather than to isolate short-term temperature effects. Consistency of temperature–mortality relationships across age groups supports the robustness of the estimated heat effects but does not by itself define their temporal scale.

Four statistical indicators were used for model evaluation: the coefficient of determination (R²), root mean square error (RMSE), relative root mean square error (rRMSE), and mean error (ME). Together, the comparative modeling framework, pre-COVID analysis, and sensitivity tests ensure that the estimated temperature–mortality relationships primarily reflect short-term temperature variability rather than long-term climatic trends, demographic shifts, or pandemic-related fluctuations. The technical workflow of the modeling and comparison scheme is summarized in Figure 2, and additional methodological details are provided in the Supplementary Materials.

Figure 2. The technical flowchart of this study.

Read detailed figure description

This flowchart details cumulative heat effects impact modeling across five steps: Data Acquisition, Preprocessing, Feature Extraction, Statistical Modeling, and Result and Analysis. Data acquisition uses MOD11A1 LST, ERA5-Land Temperature, and Eurostat mortality records. Preprocessing handles missing data via linear interpolation, validates imputation, and merges data into a CSV. Feature extraction calculates D_LST, N_LST, TDif, Tmean, and Cumulative Index PAT. Statistical modeling adjusts PAT from 100 percent to 0, determines optimal PAT, adjusts weekly data, and determines cumulative days using a random forest model with 5-fold cross validation. Results show spatial distribution at the same temperature range, accumulated temperature enhances model performance, and cumulative heat effects on mortality persist for 4 weeks. No significant difference exists between satellite and reanalysis data. A scatter plot of predicted versus true mortality shows R squared equals 0.623 and RMSE equals 32.147 deaths per million. Key findings: nighttime heat exposure and accumulated heat exposure are important for health risks, previous weeks temperature contributes 70 percent, and the model performs well on 90 independent validation sites.

Display full size

3.1 Association between temperature and mortality

3.1.1 Differences in day-night exposure response

Between 2000 and 2023, the impact of daytime and nighttime temperature exposure on mortality in Europe showed significant differences. As temperatures deviated from the MMT, the relative risk (RR) of mortality increased. The contrast between daytime and nighttime exposure revealed distinct mortality response patterns (Figure 3). During the daytime, when temperatures exceeded the MMT (31.5 °C), the RR showed a gradual, nonlinear upward trend. The increase in daytime RR was relatively moderate (95% RR = 1.76), suggesting some physiological adaptation to heat exposure under sunlight and activity cycles. In contrast, when nighttime temperatures surpassed the MMT (22.0 °C), mortality risk rose sharply (95% RR = 2.41, about 1.4 times higher than daytime). This suggests that nighttime heat has a more severe impact on mortality (Royé et al. Citation2025). While extreme nighttime heat events were rare, mortality peaked at higher temperatures, indicating that even a few extreme events contributed significantly to overall mortality (Wu et al. Citation2025). The differences in day‒night risk patterns may be linked to physiological disruptions caused by nighttime heat, such as impaired sleep quality and increased cardiovascular strain, which elevate mortality risk, especially among vulnerable populations (He et al. Citation2022). This underscores the public health risks associated with nighttime heat exposure.

Figure 3. Temperature-related risk of death during 2000–2023. (a, b) Cumulative relative risk of death (unitless) in Europe for the overall population, women, and men during daytime (a) and nighttime (b), together with their 95% CIs (shadings), and the gray curve represents the frequency distribution of temperature, the bar chart shows the number of deaths. (c, d) Cumulative relative risk of death for age 0–14, 15–59, 60–79, and 80+ during day (c) and night (d), together with their 95% CIs (shadings).

Read detailed figure description

Four-panel graph shows relative risk of death vs. temperature. Panel a: Daytime relative risk for both sexes, men, and women. Risk is highest at low temps, lowest around 30C, then increases. D_LST frequency peaks around 30C. Deaths increase from 15-40C, then decrease. Panel b: Nighttime relative risk for both sexes, men, and women. Risk is highest at low temps, lowest around 20C, then increases. N_LST frequency peaks around 20C. Deaths increase from 5-25C, then decrease. Panel c: Daytime relative risk for age groups 0-14, 15-59, 60-79, and 80+. Risk is highest at low temps, lowest around 30C, then increases. 80+ group shows highest risk at both extremes. Deaths increase from 15-40C, then decrease. Panel d: Nighttime relative risk for age groups 0-14, 15-59, 60-79, and 80+. Risk is highest at low temps, lowest around 20C, then increases. 80+ group shows highest risk at both extremes. Deaths increase from 5-25C, then decrease.

Display full size

3.2 Temperatures and heat-related mortality numbers and rates

3.2.2 Regional heterogeneity

Heat-related mortality also exhibited substantial regional heterogeneity across European cities and regions. Southern European cities generally showed higher mortality rates than northern and western Europe. For example, in Athens, the mortality rate was 318 deaths per million for women (95% CI = 256–380) and 264 deaths per million for men (95% CI = 193–336). In Rome, the corresponding rates were 359 deaths per million for women (95% CI = 284–434) and 178 deaths per million for men (95% CI = 127–229). In contrast, northern cities such as Stockholm and Oslo showed relatively lower mortality rates, with Stockholm's female mortality rate at 172 deaths per million (95% CI = 128–216) and Oslo's female mortality rate at 148 deaths per million (95% CI = 102–193). These regional contrasts suggest substantial geographic variation in heat-related health burdens.

The spatial heterogeneity was also evident in the relative risk (RR) patterns under comparable temperature conditions (Figure 4; Supplementary Material Figs. S4–S9). Even within similar daytime or nighttime temperature ranges, RR values varied considerably across regions, indicating that heat-related mortality risk is not determined solely by latitude or immediate heat intensity. For example, in Northern European coastal regions, RR was relatively high within daytime temperature ranges of 30 to 35 °C and nighttime temperature ranges above 25 °C compared with Southern European Mediterranean cities (Figs. S4–S9). In contrast, in Central and Eastern European countries such as Hungary and Romania, mortality risk rose sharply when daytime temperatures exceeded 35 °C. Southern European cities generally maintained relatively lower RR values under similar high-temperature ranges, likely reflecting greater multiweek adaptation to warm climates.

Figure 4. Relative risk (RR) for both sexes under daytime (a) and nighttime (b) across European cities during the warm season (2000–2023).

Read detailed figure description

The two-panel map shows relative risk across European cities. Panel a, labeled Daytime, displays a map of Europe with various regions outlined. A legend on the right indicates relative risk values from 1.00 to 1.10 (lightest shade) up to greater than 2.00 (darkest shade). Regions with higher relative risk are concentrated in Scandinavia, Eastern Europe, and parts of Southern Europe. Insets show magnified views of England, Poland, and Belgium and the Netherlands. Cities like Innlandet, Stockholm, Kurzeme, Klaipeda, London, Walbrzych, Koimbla, Barcelona, Rome, and Athens are labeled. Panel b, labeled Nighttime, shows a similar map of Europe with the same legend. The distribution of relative risk values appears visually similar to panel a, with some regional variations in intensity. Both panels show a general trend of higher relative risk in specific urban and surrounding areas, with the darkest shades indicating risk greater than 2.00.

Display full size

Regional heterogeneity was further reflected in differences between daytime and nighttime heat impacts. In some Northern and Eastern European cities, nighttime heat exposure was associated with a more pronounced increase in RR, whereas in Southern and Central Europe, daytime heat had a relatively stronger effect (Figure 4). The additional mortality risk associated with nighttime heat exposure across multiple regions suggests an important role of reduced nocturnal recovery and sustained physiological stress. Taken together, these findings indicate that regional differences in heat-related mortality risk are shaped not only by temperature exposure itself but also by multiweek heat exposure history, adaptive capacity, and local environmental and socioeconomic conditions across Europe.

3.3 The impact of cumulative heat effect on heat-related mortality

3.3.1 Optimal proportion of cumulative temperature

Instantaneous temperature alone cannot fully explain spatial differences in mortality rates, as prolonged high temperatures pose a continuous health threat. We hypothesize that cumulative heat effects contribute significantly to health risks, and we used machine learning models to validate this hypothesis. Models were developed based on single temperatures (D_LST or N_LST), composite temperatures (D_LST and N_LST), and cumulative heat (adding temperatures from the previous three weeks as input variables). The results show that the cumulative heat effect typically lasts three weeks. This conclusion is further supported by model performance comparisons across different cumulative exposure windows (0–3 weeks), which show consistent improvement and a plateau at approximately three weeks (Supplementary Materials Fig. S10). The model's performance is best when the proportion of cumulative temperature (PAT) reaches 70% (R² = 0.679, rRMSE = 11.21%, Figure 5a). Increasing PAT further reduces predictive accuracy, indicating that the current week's temperature contributes 70% to mortality rate predictions, with subsequent weeks weighted by exponential decay.

Figure 5. (a) Model-estimated mortality under different PAT; (b) Importance percentage of various temperature characteristics; (c) Model performance (R²) in estimating mortality across 25 typical cities with different input combinations (estimated under PAT = 70%).

Read detailed figure description

The 3-panel figure shows model performance. Panel a is a scatter plot with a fitted curve, showing R squared on the vertical axis from 0.62 to 0.69 and Proportion of accumulate temperature, PAT, on the horizontal axis from 0 percent to 100 percent. The curve rises from 0.63 at 0 percent PAT, peaks at 0.679 at 70 percent PAT, then declines to 0.66 at 100 percent PAT. Points are colored by rRMSE, ranging from 11 percent to 15 percent. Panel b is a horizontal bar graph showing Importance rate on the horizontal axis from 0 percent to 20 percent and Feature on the vertical axis. N LST has the highest importance rate at approximately 20 percent, followed by D LST at approximately 18 percent. Other features like N LST AT1, D LST AT1, LST Mean, LST Dif, N LST AT2, D LST AT2, N LST AT3, and D LST AT3 show decreasing importance rates. Panel c is a heatmap with cities on the horizontal axis and Input combination on the vertical axis, including Night3, Night2, Night1, Night0, Day3, Day2, Day1, Day0, Com3, Com2, Com1, Com0, and Average. The heatmap shows varying values from 0.2 to 0.8, with brighter regions indicating higher values.

Display full size

The day‒night composite temperature model outperforms the single-temperature model in explaining heat-related mortality (Figure 5b). Nighttime LST (N_LST) has a stronger impact than daytime LST (D_LST), emphasizing the role of nighttime high temperatures in health risks. Incorporating cumulative heat significantly improves mortality rate prediction accuracy, with the highest R² increase reaching 75% compared to models using only current-week temperatures (R² = 0.388–0.489). In 25 cities (Figure 5c), model performance varied regionally, with Milan, Amsterdam, and Bucharest performing well, while Oslo and Helsinki showed weaker results (the input combinations can be found in Supplementary Material Table S8). These findings highlight the importance of nighttime temperatures and cumulative heat in predicting health risks.

Table 2 summarizes the performance of the RF model incorporating full temperature accumulation effects. The overall accuracy for the full sample is high (R² = 0.679 ± 0.02). Model performance increases with age, reaching the highest level in the 80+ group (R² = 0.682 ± 0.02; rRMSE = 16.2% ± 1.0), while the 0–14 group shows lower accuracy (R² = 0.558 ± 0.02; rRMSE = 42.6% ± 2.0). The model performs better for women (R² = 0.687 ± 0.02) than for men (R² = 0.656 ± 0.02).

Table 2. Performance of the models—R², RMSE (deaths per million), rRMSE (%), and ME (deaths per million)—on the test (hold out) set obtained with bootstrapping using RF (PAT = 70%) (Bold in the table indicates that the statistical indicators of the model perform best in the corresponding age group or gender).

Download CSVDisplay Table

To assess whether multi-week or seasonal temporal patterns affected model performance, a detrending analysis was conducted (Supplementary Materials Table S9). The time-only model showed poor accuracy (R² = 0.214), while the temperature-only model achieved high accuracy (R² = 0.679), comparable to the temperature + time model (R² = 0.693). A shuffled temperature-only model yielded slightly lower but similar performance, supporting the stability of the temperature–mortality relationship. Additionally, a pre-COVID model trained using data from 2000–2019 produced nearly identical performance to the full-sample model, suggesting that pandemic years did not bias the results. The age-stratified pattern of increasing accuracy further suggests that the stronger performance in older groups is consistent with their higher observed heat vulnerability, but does not by itself establish a specific physiological mechanism or exclude other age-related factors. Overall, these analyzes confirm that model performance is primarily driven by temperature predictors, and that temporal trends, demographic shifts, and pandemic-related mortality fluctuations have limited influence on the estimated temperature–mortality relationship.

3.3.2 Spatial distribution of model performance improvement

Nighttime LST exerts a significantly greater impact on model accuracy than daytime (Supplementary Materials Fig. S11), indicating that nighttime LST poses a more severe threat to thermal health. When both daytime and nighttime LST are considered simultaneously, the model temperatures increase further (Figure 6a). At the aggregate level, the baseline model without cumulative temperature variables shows moderate performance across the 122 cities, with R²generally in the range of 0.40–0.50 and minimum rRMSE values around 13%–15%, indicating that contemporaneous temperature alone explains only part of the mortality variability. Nevertheless, models relying solely on weekly temperatures show limited accuracy and clear regional variability. Northern Europe (e.g. Sweden, Norway), Western Europe (e.g. Belgium, Netherlands), and Southern European coastal cities (e.g. Spain, Portugal) generally exhibit higher R² values (often exceeding 0.60 or 0.80), indicating sensitivity to short-term temperature fluctuations. In contrast, some inland and Central Eastern European cities (e.g. Poland, Austria, and the Czech Republic) display lower R²(below 0.30), suggesting that single-week temperatures fail to capture cumulative heat impacts. The spatial distribution of rRMSE (Figure 6b) further indicates that low R² does not necessarily imply large prediction errors. In several cities, a relatively low R² is accompanied by modest rRMSE, likely due to limited variability in the mortality series. Therefore, model performance should be evaluated jointly using both R² and rRMSE. Figure 6c illustrates the contribution of different temperature variables. Including both D_LST and N_LST improves R² and reduces RMSE, highlighting the importance of incorporating both daytime and nighttime temperatures. The corresponding rRMSE patterns support this result, showing reduced prediction errors across many regions, even where R² gains are modest. Nighttime temperature contributes slightly more to model accuracy, reinforcing its importance in capturing heat-related effects.

Figure 6. Performance of the trained RF models (the cumulative temperature variable was not considered). Spatial distribution of R² (a) and rRMSE (b) of the model in different cities using both D_LST and N_LST as input. (c) R² and RMSE of the model with the above three inputs.

Read detailed figure description

The three panel figure shows R squared and rRMSE of mortality estimates across Europe and performance metrics for different inputs. Panel a displays a map of Europe with R squared values for mortality estimates. A legend indicates R squared values from 0.10 to 0.20, 0.20 to 0.30, 0.30 to 0.40, 0.40 to 0.50, 0.50 to 0.60, and greater than 0.60. Insets show magnified views of England, Poland, and Belgium and the Netherlands. Panel b displays a map of Europe with rRMSE percentage values for mortality estimates. A legend indicates rRMSE values from 0.0 to 7.0, 7.0 to 9.0, 9.0 to 11, 11 to 13, 13 to 15, and greater than 15. Insets show magnified views of England, Poland, and Belgium and the Netherlands. Panel c contains two grouped scatter plots. The top plot shows R squared values on the vertical axis from 0.0 to 0.6. The horizontal axis shows three groups: dayLST, nightLST, and dayLST plus nightLST. Each group has subcategories of Total, Man, and Woman. The R squared values generally range from 0.3 to 0.5 across all categories. The bottom plot shows rRMSE percentage values on the vertical axis from 0 to 30. The horizontal axis has the same three groups and subcategories. The rRMSE values generally range from 10 to 20 percent across all categories.

Display full size

When cumulative temperature variables are included (Figure 7a), model performance improves markedly across most regions. In the full model, R² increases to approximately 0.65–0.68, while the minimum rRMSE decreases to about 10%–11%, demonstrating a clear enhancement compared with the baseline model. The gains are most evident in Northern, Central, and some Southern European coastal regions, where R² increases substantially, with some cities exceeding 0.80. This pattern reflects the importance of multiweek heat exposure, particularly in regions with lower average temperatures and fewer extreme events. In contrast, improvements are smaller in the UK, parts of Southern Europe (e.g. Spain, Italy), and Eastern Europe (e.g. Romania), likely due to climatic conditions and adaptive capacity. The spatial distribution of rRMSE (Figure 7b) shows an overall reduction in prediction error after incorporating cumulative variables. From a gender perspective (Supplementary material Fig. S12), cumulative heat effects further improve prediction accuracy for female mortality, particularly in Northern and Central Europe, where rRMSE decreases (Figure 7c). Overall, cumulative temperature enhances model performance across most regions, while more limited improvements in some areas likely reflect climate adaptation and lower variability. These results demonstrate that cumulative heat exposure improves both explanatory power and prediction stability, highlighting its importance for heat–health risk assessment and supporting the need for region-specific adaptation strategies.

Figure 7. Performance of the trained RF models considering temperature accumulation effects. Spatial distribution of R² (a) and rRMSE (b) for models estimating mortality rates for the total population (both sexes). (c) R² and RMSE of the model for estimating mortality rates for the total population, man, and woman.

Read detailed figure description

The two-panel figure shows a map and two scatter plots. Panel a displays a map of Europe with regions shaded by R squared values. The legend indicates R squared values less than 0.40 are light blue, 0.40 to 0.50 are blue, 0.50 to 0.60 are light red, 0.60 to 0.70 are red, 0.70 to 0.80 are dark red, and greater than 0.80 are very dark red. Insets show magnified views for England, Poland, and Belgium and the Netherlands. Panel b contains two scatter plots. The left scatter plot shows R squared on the vertical axis from 0.0 to 1.0, with data points for Total, Man, and Woman categories on the horizontal axis. For all three categories, R squared values range from approximately 0.4 to 0.8, with mean values around 0.6 to 0.7. The right scatter plot shows rRMSE percentage on the vertical axis from 8 to 20, with data points for Total, Man, and Woman categories on the horizontal axis. For all three categories, rRMSE percentage values range from approximately 8 to 20, with mean values around 12 to 14.

Display full size

In summary, the baseline model without cumulative temperature variables shows moderate predictive performance, with R²generally ranging from 0.40 to 0.50 and minimum rRMSE values around 13%–15%, while performance for men is consistently slightly lower than for women, indicating higher predictive uncertainty in male mortality estimates. After incorporating cumulative temperature variables, R² rises to approximately 0.65–0.68, and minimum rRMSE decreases to about 10%–11%, corresponding to an average improvement of 36%–63% in R² and a reduction of 23%–26% in rRMSE across the 122 cities. These systematic gains, observed across regions and demographic groups, provide strong empirical support for the cumulative heat framework. They demonstrate that multiweek cumulative heat exposure substantially enhances model explanatory power and prediction accuracy, indicating that heat-related mortality is influenced not only by short-term temperature fluctuations but also by the persistence and accumulation of heat stress over time. This highlights the necessity of explicitly accounting for cumulative exposure in large-scale heat–health assessments.

Finally, to assess potential uncertainties arising from cloud-induced data gaps in satellite-derived MODIS LST, we conducted a parallel analysis using seamless temperature data from the ERA5-Land reanalysis. The RF model achieved an R² of 0.623 and an RMSE of 32.147 deaths per million when driven by ERA5-Land temperatures, demonstrating a high degree of consistency with the results obtained using MODIS LST data (Supplementary Material Fig. S13). This strong agreement between two independent temperature sources provides a robust validation of the MOD11A1 dataset's reliability for estimating heat-related mortality. And it also confirms that the identified cumulative heat effects are a robust phenomenon, not an artifact of the specific temperature data source used.

3.4 Differences in cumulative heat effects across regions and independent verification

At the city scale, the duration of cumulative heat effects varies across regions (Figure 8a). Cumulative heat effects are longer in Northern and Southern Europe, with heat exposure lasting four to five weeks in the north and slightly shorter durations in the south. In contrast, Central and Eastern European inland areas (e.g. Poland, Czech Republic, and Austria) experience shorter cumulative effects, typically lasting two to three weeks. The longer duration in the northern and southern regions is largely due to extended periods of extreme heat events, while Central Europe experiences more fluctuating, shorter heat spells. Regions in the north and south, with more extreme seasonal temperature variations, face prolonged heat exposure, while inland areas in Central Europe are exposed to more transient heat events, limiting the persistence of cumulative effects. Additionally, more developed economic development can also help to prolong treatment time. Despite shorter durations, frequent heat exposure still poses significant health risks.

Figure 8. (a) Spatial Distribution of cumulative heat exposure weeks in European cities; (b) estimated mortality rate (deaths per million) by replacing 2019 temperature data with 2022 temperature data. Weekly high-temperature data from 2022 was replaced with the corresponding weeks from 2019, and a machine learning model was used to predict mortality based on the modified data (with temperature data from the first three weeks of 2019 and the corresponding weeks from 2022 as inputs). The optimal cumulative heat exposure duration for each city was determined using a stepwise sensitivity analysis, in which lagged temperature windows were progressively extended, and the window yielding the best model performance (R² and rRMSE) was selected.

Read detailed figure description

The two-panel figure shows a map and a box plot. Panel a is a map of Europe displaying cumulative heat exposure weeks in various cities. A legend indicates shades from 1 Week to 5 Weeks. Insets show magnified views for England, Poland, and Belgium and the Netherlands. Panel b is a box plot with a vertical axis labeled Deaths per million, ranging from 0 to 400. The horizontal axis shows three categories: 2019, 2022, and 2022 replace 2019. The 2019 box plot has a median of 175.8. The 2022 box plot has a median of 194.2. The 2022 replace 2019 box plot has a median of 184.7. The median for 2022 is higher than 2022 replace 2019, which is higher than 2019.

Display full size

To assess the impact of temperature changes on mortality, we substituted the 2022 weekly high-temperature data with 2019 data. The predicted mortality rate after the substitution was 184.7 per million, between the 2019 rate of 175.8 and the 2022 rate of 194.2 per million (Figure 8b), indicating the sensitivity of mortality to interannual differences in temperature conditions, particularly under years with more frequent or intense heat events. However, this comparison primarily reflects differences in overall temperature magnitude between the two years and does not, by itself, isolate the effect of cumulative heat exposure duration. To further examine this mechanism, we conducted a simulation analysis using four synthetic exposure scenarios (A–D; Supplementary Material Table S10). These scenarios do not represent four separate models, but rather four designed temporal patterns of heat exposure that were input into the same trained RF model. Scenario A (continuous heat exposure) resulted in higher mortality than Scenario B under comparable heat conditions, indicating that more persistent heat exposure increases mortality risk. Scenarios C and D, which represent alternating or more sporadic heat events, showed that intervals between heat events weakened the cumulative heat effect and reduced the average mortality during heat-affected weeks. These controlled simulations provide more direct evidence that the temporal persistence and continuity of heat exposure contribute to cumulative heat effects, independent of overall temperature magnitude. Therefore, while the year-substitution analysis illustrates the overall influence of heat conditions, the role of cumulative heat duration is more robustly supported by the synthetic scenario experiments.

For model validation, we selected 90 cities across Europe. The RF model performed well with new data, with performance similar to the original cities (Figure 9a–c). However, gender differences were noted. For men, R² values exceeded 0.6 only in northern Norway and a few southern cities, while most regions had values around 0.5. The model's accuracy for women was significantly better, with high performance in southern cities, such as in Spain and Italy, and lower R² values only in the UK and inland cities. The final model achieved an R² of 0.584 in the independent cities, a 14% decrease, with an RMSE of 35.276, a 31.5% increase (Figure 9d).

Figure 9. Performance of the trained RF models considering temperature accumulation effects at 90 independent validation cities (using MOD11A1 temperature data). Spatial distribution of R² for models estimating mortality rates for the total population (a), men (b), and women (c) separately. (d) Scatter plot of predicted mortality and actual mortality using the RF model at 90 independent validation cities.

Read detailed figure description

The four panel figure, arranged two by two, shows R squared spatial distribution for mortality rates and a scatter plot. Panel a, titled Estimation of total mortality rate, displays a map of Europe with regions shaded based on R squared values, ranging from less than 0.40 to greater than 0.80. Inset maps show England and Belgium and Netherlands. Panel b, titled Estimation of men mortality rate, shows a similar map with R squared values for men. Panel c, titled Estimation of women mortality rate, presents another map with R squared values for women. Panel d, titled Performance of estimated mortality rate, is a scatter plot with True mortality on the horizontal axis and Predicted mortality on the vertical axis, both ranging from 0 to 800. The plot shows a dense cluster of points along a diagonal line, indicating a strong correlation. A solid diagonal line represents the ideal prediction, and a dashed line shows the actual trend. The R squared value is 0.584 and RMSE is 35.276 death per million.

Display full size

4.1 Modeling approaches and comparative performance

When assessing the impact of heat on mortality, model selection is crucial, including handling nonlinear relationships, considering cumulative heat effects, and ensuring regional comparability (Cui et al. Citation2023; Huang et al. Citation2025). Traditional Poisson regression models are widely used due to their simplicity, but they assume linear relationships, making it difficult to capture complex nonlinear effects and interactions (Weaver et al. Citation2015). Additionally, Poisson regression struggles with collinearity, cumulative effects, and overfitting, limiting predictive accuracy (Lüthi, Fairless, and Fischer Citation2023). In contrast, machine learning models can better capture nonlinear relationships, handle complex interactions, and improve prediction (Yin et al. Citation2023). For instance, in this study, RF demonstrated the best balance between predictive performance, robustness, and interpretability, outperforming penalized regression (PR) and XGBoost (Supplementary Material Table S11). While PR improves stability over traditional regression models by addressing multicollinearity through regularization, it still struggles to model complex nonlinear interactions (Bellavia Citation2025; Kim and Kim Citation2022). On the other hand, XGBoost showed strong predictive accuracy (Bora et al. Citation2025), but its greater sensitivity to hyperparameter settings may increase the complexity of model tuning. Therefore, in our heat–mortality analysis, we selected random forest (RF) because, within our modeling framework, it allowed more straightforward interpretation of feature contributions and cumulative heat patterns, and facilitated more consistent comparison across cities and population groups. RF, however, consistently delivered robust performance across regions and maintained interpretability through permutation-based feature importance, making it particularly suitable for large-scale, comparative heat–mortality analysis (Shao, Ahmad, and Javed Citation2024). Its ability to model both the short-term and cumulative effects of heat exposure across diverse regions supports its selection as the primary model in this study.

Most existing studies use city-specific models (Anderson and Bell Citation2009; Chestnut et al. Citation1998; Kephart et al. Citation2022; Wang et al. Citation2024), which, while reflecting local characteristics, make cross-city comparisons difficult and fail to quantify cumulative heat effects (Pattenden, Nikiforov, and Armstrong Citation2003). Our results show that using a unified modeling framework across all cities improves comparability and significantly boosts model accuracy (Supplementary Material Fig. S14), with an R² of 0.679, a 95.1% improvement over traditional Poisson regression (R² = 0.348), and a reduced rRMSE of 11.21%. These improvements enhance confidence in the robustness of the modeling framework and support its suitability for assessing cumulative temperature–mortality associations across regions. Overall, the DLNM-based analysis and RF-based modeling provide consistent evidence that heat exposure is strongly associated with increased mortality risk, while emphasizing different dimensions of the relationship. The epidemiological model quantifies interpretable relative risk patterns across temperature ranges and demographic groups, whereas the ML model captures the temporal persistence of heat impacts by explicitly incorporating lagged accumulated temperatures. Importantly, regions showing higher relative risks under nighttime heat conditions also tended to exhibit longer cumulative heat durations and stronger contributions from prior-week exposure, supporting the interpretation that nighttime heat acts as a key driver of cumulative thermal stress.

4.2 Cumulative and nighttime heat effects

Previous studies have shown that mortality risks associated with nighttime high temperatures are significantly higher than those linked to daytime heat (Liu et al. Citation2025; Wu et al. Citation2023). Elevated nighttime heat is strongly associated with disrupted sleep and circadian temperature regulation (Obradovich et al. Citation2017), impaired immune responses (Bryant, Trinder, and Curtis Citation2004; Meier-Ewert et al. Citation2004), and potential declines in psychological and cognitive performance (Waters and Bucks Citation2011). Consistently, we found that the relative risk of mortality under nighttime heat exposure (RR = 1.62, 95% CI = 1.45–1.78) was substantially higher than that under daytime high temperatures (RR = 1.37, 95% CI = 1.23–1.52), confirming a strong statistical association between nighttime heat and elevated mortality (Royé et al. Citation2025). This disparity may arise because daytime heat can often be mitigated through reduced outdoor activity or cooling measures, whereas nighttime heat is more difficult to avoid (He et al. Citation2022; Rippstein et al. Citation2023). Nighttime temperature exposure, therefore, reflects living environment quality and socioeconomic constraints (Romitti et al. Citation2025), as economically underdeveloped areas frequently lack adequate cooling facilities. For example, air-conditioning adoption and usage remain low in many European countries, even in relatively hot regions such as Spain and Italy (Randazzo, De Cian, and Mistry Citation2020). Sustained nighttime heat also reduces diurnal temperature range, intensifies sleep disturbance and cumulative physiological strain, and may increase vulnerability among sensitive groups. Notably, older adults (>60 years) are particularly susceptible, with relative risks reaching 2.25 (95% CI = 1.98–2.57).

Focusing solely on daytime or nighttime temperatures is insufficient to capture the multiweek health impacts of high-temperature exposure. Temperature-related mortality involves at least two distinct temporal components: short-term lag effects (Gasparrini Citation2014) and multiweek cumulative effects (Liu et al. Citation2025). Lag effects, widely examined using Distributed Lag Nonlinear Models (DLNM), describe delayed mortality responses occurring within several days after extreme events and reflect immediate but time-shifted physiological stress (Gasparrini, Armstrong, and Kenward Citation2010). In our framework, the cumulative effects quantified in this study represent the aggregated influence of sustained high-temperature exposure over multiple consecutive weeks, and may be interpreted as an extension of the temporal exposure structure beyond short-term lag responses rather than a completely separate mechanism. Our results indicate that cumulative effects of prolonged heat exposure are a more critical driver of heat-related mortality, reflecting the potential compounding physiological burden and the gradual weakening of thermoregulatory capacity under persistent environmental strain (Shin et al. Citation2024). We found that high-temperature exposure during the first three weeks significantly increased mortality risk in the current week, with exposure from the previous week contributing up to 70% of the current-week risk. Although the impact declines exponentially over time, it remains significant for an average of four weeks, far longer than the 3–7 day lag reported in previous studies (Wu et al. Citation2024). These findings highlight that lag and cumulative effects capture different aspects of temporal dependence in heat-related mortality: lag effects emphasize immediate delayed responses, whereas cumulative effects capture the persistence and accumulation of thermal stress over extended periods. Incorporating cumulative temperature variables also substantially improved predictive performance, with an average R² increase of 15%, further confirming the added value of considering cumulative heat effects. This cumulative dimension provides new insights into multiweek heat stress beyond short-term lag patterns, particularly in economically and medically resource-limited regions where prolonged exposure may amplify health risks.

This study found that the duration of cumulative heat effects varies markedly across Europe. Both Southern Europe (e.g. Spain, Portugal) and Northern Europe (e.g. Norway, Sweden) exhibit longer-lasting cumulative heat effects, with impacts in many cities persisting for four weeks or more. In contrast, Central inland regions (e.g. Germany, Poland, Czech Republic) show shorter durations, typically one to two weeks. Several factors may explain this pattern. Southern Europe's consistently high temperatures and frequent heatwaves create sustained thermal stress, while Northern Europe's relatively cool climate means populations are less physiologically adapted to heat; when extreme events occur, the effects linger longer before recovery (Ballester et al. Citation2023). Central Europe, with more moderate and fluctuating summer temperatures, experiences shorter heat episodes, resulting in weaker cumulative effects. Moreover, coastal regions in the north and south tend to retain heat longer due to higher humidity and slower nighttime cooling, reinforcing persistence in cumulative impacts. In contrast, inland regions cool more rapidly after sunset, interrupting continuous heat accumulation. Urban characteristics, such as dense built environments and limited vegetation in southern and northern cities, further amplify nighttime heat retention (Estoque et al. Citation2020; Relvas et al. Citation2025), intensifying long-duration health effects.

4.3 Regional heterogeneity and adaptive capacity

Socioeconomic conditions and infrastructure substantially influence the persistence of heat effects (Iungman et al. Citation2023). In wealthier Western European cities, well-developed cooling systems, healthcare services, and urban planning can mitigate prolonged heat risks and shorten the duration of cumulative effects (Hone et al. Citation2020). In contrast, economically underdeveloped regions may experience multiweek-lasting impacts, placing heavier burdens on healthcare systems and extending health threats (Heidenreich and Thieken Citation2024). To further characterize these regional disparities, we quantified the distribution and variability of key socioeconomic and educational indicators across European regions (Supplementary Material Fig. S15), including per capita GDP, medical GDP, and tertiary education levels (available at https://ec.europa.eu/). The high coefficients of variation indicate substantial regional heterogeneity in economic resources and adaptive capacity, which may help explain the observed differences in cumulative heat exposure duration.

Adaptive capacity modifies cumulative heat vulnerability through several interrelated mechanisms. First, access to cooling infrastructure (e.g. air conditioning, building insulation, and urban green spaces) can reduce both the intensity and persistence of individual heat exposures, thereby limiting the accumulation of physiological stress across consecutive days or weeks (Basu and Samet Citation2002; Hondula et al. Citation2015). Second, healthcare accessibility and early warning systems enable more timely intervention during prolonged heat events, reducing delayed or cumulative mortality risks (Ebi et al. Citation2021). Third, socioeconomic conditions influence behavioral adaptation, such as the ability to modify activity patterns or seek cooler environments, which further mitigates sustained exposure (Sera et al. Citation2019). Together, these factors determine how quickly populations can recover from heat stress and thus directly influence the duration and magnitude of cumulative heat effects. Population adaptive capacity further contributes to these differences: communities exposed to high temperatures over long periods may develop stronger physiological adaptation, whereas high-latitude populations experiencing fewer extreme events remain more vulnerable. This interaction between environmental exposure history and socioeconomic capacity suggests that cumulative heat effects are not solely driven by temperature magnitude, but also by the ability of populations to buffer, respond to, and recover from repeated heat stress over time. Thus, cumulative-effect duration can vary even under similar heat events, underscoring the need for region-specific heat protection measures.

4.4 Methodological implications of multiweek cumulative heat exposure

In this study, we applied an RF model together with the concept of multiweek cumulative heat exposure to address an important gap in understanding how sustained heat may shape mortality risk beyond conventional short-lag structures. Existing temperature–mortality studies, particularly those based on time series or DLNM frameworks, have provided robust evidence on nonlinear and delayed associations, but they have predominantly been developed at the daily scale and are typically interpreted within lag windows of several days to a few weeks (Gasparrini et al. Citation2015). Similarly, heatwave studies have shown that heat intensity and event duration are important for mortality risk, but most of this literature operationalizes duration as consecutive hot days rather than as cumulative exposure persisting across multiple weeks (Guo et al. Citation2017; Xu et al. Citation2016). Accordingly, our use of the term cumulative heat effect refers specifically to a medium-term temporal dimension—namely, the persistence and carry-over of heat-related mortality risk across consecutive weeks—rather than to chronic health burdens at seasonal, annual, or interannual scales. This distinction is important because sustained heat exposure may not only trigger acute responses after individual hot days, but may also progressively increase vulnerability when recovery between hot periods is incomplete. In this context, nighttime heat may be especially relevant because elevated nocturnal temperatures reduce the opportunity for physiological cooling and recovery, thereby prolonging internal heat strain and potentially amplifying mortality risk over subsequent weeks (He et al. Citation2022). By incorporating both heat intensity and temporal persistence, our cumulative heat framework extends the interpretation of heat-related risk from short-term lagged responses to a broader medium-term accumulation process. As an ensemble learning method, RF is well suited to this purpose because it can flexibly model nonlinear, high-dimensional, and regionally heterogeneous heat–mortality relationships without imposing strong distributional assumptions (Breiman Citation2001). Its ensemble structure extracts long-term signals from smoothed weekly series, while its nonparametric nature reduces bias from temporal aggregation (Ballester et al. Citation2024). Combined with cross-validation and feature selection, RF enables robust quantification of cumulative heat impacts across regions and populations, maintaining high predictive accuracy while avoiding overfitting (Aldrich Citation2020). Therefore, rather than suggesting that conventional exposure–mortality models are inadequate, our intention is to show that the present framework complements existing epidemiological approaches by capturing an additional temporal dimension of heat-related mortality that is less directly represented in standard daily lag models.

From a methodological perspective, DLNM remains the most widely adopted epidemiological framework for estimating nonlinear and lagged temperature–mortality associations, offering clear interpretability in terms of exposure–response curves and relative risks (Gasparrini, Armstrong, and Kenward Citation2010). However, DLNM-based analyzes are typically constrained by predefined lag structures and may be less flexible when modeling heterogeneous cumulative effects across multiple regions. In contrast, ML models such as RF can integrate multiple temperature indicators (e.g. daytime/nighttime extremes, diurnal range, and accumulated anomalies) and capture complex nonlinear interactions without strong parametric assumptions. By combining these two approaches, this study provides both interpretable epidemiological evidence and a data-driven quantification of cumulative heat persistence, which improves the robustness of conclusions regarding prolonged heat exposure.

The cumulative heat effect concept offers a new perspective on multiweek heat-related mortality, particularly under global warming and urbanization, where cities increasingly experience more frequent and prolonged heat exposure beyond short heatwaves, spanning multiple seasons or years (Lee et al. Citation2020). This framework more accurately captures population health burdens, especially among vulnerable groups such as the elderly, as prolonged heat stress can gradually weaken physiological regulatory mechanisms and increase disease and mortality risks (Bi et al. Citation2023). By integrating the RF model with cumulative heat effects, this study systematically quantifies multiweek mortality risks across regions and demographic groups. Importantly, the novelty of this study does not lie simply in replacing ERA5-Land with MODIS LST, nor in claiming that LST is physically equivalent to air temperature. Rather, the key contribution lies in using satellite-derived daytime and nighttime LST as a high-resolution, spatially explicit exposure framework to quantify cumulative heat effects across a large and climatically heterogeneous set of European cities, while using ERA5-Land as an independent reference dataset to test the robustness of the main findings. Although satellite-derived land surface temperature (LST) differs physically from near-surface air temperature experienced by humans, the broadly consistent results obtained from MODIS LST and ERA5-Land indicate that the identified cumulative heat–mortality relationships are not dependent on a single temperature product. The slightly better performance of the MODIS-based model should therefore be interpreted cautiously; its added value lies primarily in its finer spatial resolution and stronger ability to capture spatial thermal heterogeneity, rather than in a dramatic gain in predictive accuracy alone. This validation demonstrates that satellite-based assessment can complement conventional reanalysis products in epidemiological research, particularly for intercity comparison and for identifying spatially heterogeneous urban heat exposure patterns that may be smoothed in coarser-resolution datasets. Moving beyond traditional short-term exposure frameworks, this approach provides a stronger foundation for public health strategies addressing the urban heat island (UHI) effect and related health threats (Watkins, Palmer, and Kolokotroni Citation2007). Overall, combining RF modeling with cumulative heat assessment offers valuable insights into heat–mortality relationships and informs effective urban health interventions.

4.5 Contextualization and contributions

To contextualize our results, substantial evidence from previous studies also links prolonged high temperatures to increased mortality. A multicountry analysis across 18 nations found that heatwaves lasting 3–4 days significantly elevated mortality risk (Guo et al. Citation2017). In Europe, the 2022 five-week summer heatwave caused an estimated 61,672 excess deaths, with mortality patterns closely aligning with cumulative thermal deviations (Ballester et al. Citation2023). Zhao et al. (2024) reported an annual average of 2.36 heatwave-related deaths per million people from 1990 to 2019, indicating a persistent long-term burden. In Germany, high-resolution spatial modeling estimated nearly 48,000 heat-related deaths between 2014 and 2023, including over 1,100 deaths during a single week-long heatwave in July 2023 (Wang et al. Citation2024). Similarly, multi-day heatwaves significantly increased mortality in Switzerland, particularly among individuals aged 80 and above (Konstantinoudis et al. Citation2018), and the 1995 Chicago heatwave showed that even five days of extreme heat could cause several hundred deaths (Semenza et al., Citation1996).

This study builds on previous research, introducing key innovations. First, we define cumulative heat exposure using a weekly framework, summing deviations from the multiweek mean temperature over the past three weeks. This medium-term metric better reflects the persistent effects of heat stress than the daily metrics used in previous studies (Vicedo-Cabrera et al. Citation2018). Second, by applying high spatial resolution (NUTS-3) and standardizing mortality data by age and sex (weekly deaths per million), we enhance the precision of inter-regional and demographic comparisons. Third, we integrated multiple data sources, including Eurostat, the UK Office for National Statistics (ONS), and national statistical offices. This approach aligns with recent advances in environmental epidemiology (Robben et al. 2025) and improves the reliability of our estimates. Our results support and expand upon global assessments by the WHO, which reported an average of 489,000 heat-related deaths annually from 2000 to 2019, with Europe accounting for 36% of the global burden (WHO, 2021). Murage et al. (Citation2024) emphasized that aging populations, such as in the UK, face disproportionately higher risks under future climate scenarios. In the US, analyzes found that homeless populations experience heat-related emergency visits at rates 27 times higher than the general population (Weckstein et al., Citation2025), highlighting the need for more granular, population-specific risk assessments. Our study, which includes age- and sex-stratified analyzes, identifies vulnerable subpopulations and quantifies their differential heat sensitivity across regions.

Substantial progress has been made in characterizing heat–mortality relationships, particularly through distributed lag nonlinear models (DLNM), which capture nonlinear and delayed effects beyond traditional linear approaches (Gasparrini Citation2021). Meta-regression methods have further synthesized city-specific risks and demonstrated how urban characteristics shape heat-related mortality, facilitating intercity comparisons (Gasparrini and Armstrong Citation2013; Schinasi et al. Citation2023; Song et al. Citation2024). However, most existing work remains localized, and large-scale harmonized transnational assessments of cumulative heat effects are still limited. This study addresses this gap by combining extensive mortality records from multiple European cities with satellite-derived land surface temperature (LST) data and RF modeling to quantify cumulative heat impacts and their spatial variability across urban settings. This framework enables broader intercity comparisons and provides new insights into interactions between sustained heat exposure and urban heat island (UHI) effects. Importantly, few studies have examined the cumulative burden of prolonged and nighttime heat; by distinguishing between daytime and nighttime conditions, we offer empirical evidence of their differing health risks. Overall, this large-scale, data-driven approach advances understanding of heat-related mortality under heatwaves and persistent warming, providing a scientific basis for anticipating and mitigating climate change–driven public health impacts and informing adaptive heat-response policies globally.

4.6 Limitations and future directions

However, several limitations remain. (1) The weekly temporal resolution of mortality data restricts the model's ability to capture short latent effects such as acute heat stress, potentially underestimating impacts from extreme heat events lasting only hours to three days. Weekly aggregation may smooth short-term variability and obscure immediate mortality responses. (2) Although mortality coverage is nearly complete, minor inconsistencies across national statistical agencies and limited imputation of missing values may introduce small biases. (3) The 1 km spatial resolution of MODIS LST limits detection of intraurban thermal heterogeneity (e.g. contrasts among built-up areas, green spaces, and waterfront zones), potentially underestimating fine-scale microclimatic effects. (4) Because the dataset emphasizes large metropolitan regions, an urban bias may exist, possibly overestimating heat-related mortality due to both urban heat island influences and the concentration of medical services that attract more recorded deaths. (5) As the model is calibrated on historical data, it may underestimate future risks from increasingly extreme and prolonged heat events; moreover, the effects of population aging and heightened vulnerability among older adults are not fully accounted for. (6) Finally, while MODIS-derived LST captures spatial thermal heterogeneity, it represents land surface conditions rather than the 2 m air temperature, apparent temperature, or indoor thermal exposure directly experienced by individuals. Therefore, LST should be interpreted as a spatial proxy for ambient heat exposure, and its epidemiological value in this study is strengthened by the consistency of results obtained from the parallel ERA5-Land analysis. 7) In addition, although machine learning models such as random forests offer strong predictive performance and flexibility in capturing nonlinear and cumulative relationships, they also have inherent methodological limitations. In particular, RF is primarily a data-driven predictive tool and does not explicitly model causal mechanisms; therefore, the identified associations between cumulative heat exposure and mortality should not be interpreted as direct causal effects. Moreover, despite the use of feature importance metrics, the internal decision-making process of ensemble models remains relatively opaque compared to traditional statistical models, which limits their interpretability and may pose challenges for directly translating model outputs into mechanistic explanations or policy-relevant causal insights.

Future research can deepen in the following directions: 1) Integrating observations from meteorological stations, wearable device data, and other multisource data (such as building insulation performance, accessibility to green spaces, and healthcare coverage) to create a socioecological exposure model, exploring climate adaptation differences. 2) Leveraging Earth Observation Systems (EOS) to enhance heat-health research, particularly in tracking real-time environmental conditions and providing comprehensive, global-scale data on heat exposure. This can improve our understanding of heat impacts across diverse geographies and inform better-targeted interventions. 3) Developing models that couple temperature and socioeconomic indicators to quantify the cost‒benefit ratios of heat mitigation measures and identify priority intervention areas. In addition, future studies should explicitly examine the relationship between socioeconomic development and the duration of cumulative heat exposure, as well as evaluate the potential moderating role of socioeconomic factors in shaping heat-related mortality risks. Such analyzes would help clarify the mechanisms underlying regional heterogeneity in cumulative heat effects. 4) Expanding research to other regions to test the cumulative effects model's generalizability and investigate how humid heat interactions amplify health risks. 5) Moreover, future studies could explore ways to enhance the interpretability of machine learning models, such as employing explainable AI (XAI) techniques, which would allow for a more transparent understanding of how different temperature variables contribute to heat-related mortality risk. These improvements will help shift heat health risk assessments from relying on a single meteorological indicator to a dynamic, multidimensional early warning system, providing more accurate scientific support for climate adaptation policies.