the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 11 Aug 2026
Long-Term Satellite Analysis of Summer Heatwaves and Heat Vulnerability in Romania (2000–2025)
Abstract. Heatwaves represent one of the most significant climate-related hazards affecting human health, ecosystems, and socioeconomic activities. This study presents the first multi-decadal, national-scale assessment of summer heatwaves and heat vulnerability in Romania using a 26-year archive (2000–2025) of MODIS Land Surface Temperature (LST) integrated with official meteorological warnings and high-resolution population-grid data. Daily satellite observations were processed to analyze long-term trends, anomalies, and consecutive hot days. A 1-km Heat Vulnerability Index (HVI) was developed by fusing satellite-derived thermal hazard data, cumulative heat-warning severity, population density, and demographic sensitivity.
The results indicate a severe and widespread intensification of summer thermal conditions across Romania, with August exhibiting the strongest surface warming trend (1.76 °C decade⁻¹), followed by July (0.82 °C decade⁻¹), and a seasonal average increase of 0.84 °C decade⁻¹. Maximum LST frequently exceeded 50 °C in the southern lowlands during the most extreme summers on record (2007, 2012, 2022, 2024, and 2025), with 2024 registering a record-breaking 56 heatwave days. The southern lowlands and the Bucharest metropolitan area emerged as the principal national hotspots of heat-related risk. Nationally, over 3.15 million inhabitants (16.5 % of the total population) reside in areas classified under high and very high vulnerability classes, with Bucharest alone concentrating 1.82 million vulnerable citizens. The proposed framework provides an operational, high-utility tool for geographic screening, directly supporting evidence-based public health interventions and climate adaptation strategies.
- Preprint
(7681 KB) - Metadata XML
-
Supplement
(6444 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-4415', Anonymous Referee #1, 01 Sep 2026
-
AC2: 'Reply on RC1', Anisoara Irimescu, 16 Sep 2026
We sincerely thank the Reviewer for the careful evaluation of our manuscript and for the constructive comments and suggestions. We have carefully considered all the points raised and revised the manuscript accordingly.
Major Concerns
1. Inhomogeneous satellite time series (Terra vs. Aqua overpass times) and 2002 data gap
Response: We thank the Reviewer for raising this important point concerning the temporal consistency of the MODIS LST record. We agree that Terra and Aqua have different nominal local overpass times and that their absolute daytime LST observations should therefore not be assumed to be directly interchangeable without considering the diurnal temperature cycle.
We have revised the manuscript to clarify the exact temporal composition of the MODIS dataset. Terra MOD11A1 observations were used for 2000–2001 and June 2002, while Aqua MYD11A1 observations were used from July 2002 onwards. Therefore, June 2002 was not omitted from the summer analysis, as might have been inferred from the previous description of the dataset. The corresponding description in Section 2.1 has been corrected accordingly.
In response to the Reviewer’s concern regarding the temporal inhomogeneity introduced by the Terra–Aqua transition, we substantially revised the trend analysis. The previously reported OLS trends based on the combined record were replaced by a non-parametric analysis based exclusively on Aqua MODIS observations. Sen’s slope was calculated on a pixel-by-pixel basis, and statistical significance was independently assessed using the Mann–Kendall test (p < 0.05). July and August trends were calculated for 2002–2025, while June and the complete seasonal JJA analysis were calculated for 2003–2025, thereby ensuring that the trend assessment is based on a temporally consistent Aqua record and is not affected by the Terra–Aqua transition. The corresponding Methods, Results, Abstract, Conclusions, and Fig. 5 have been revised accordingly.
The revised Aqua-only analysis yields mean Sen’s slopes of 0.17 °C decade⁻¹ in June, 0.83 °C decade⁻¹ in July, 2.26 °C decade⁻¹ in August, and 1.10 °C decade⁻¹ for JJA. Positive trends occur over 65%, 88%, and 100% of Romania in June, July, and August, respectively, and over 98% for JJA. The Mann–Kendall analysis indicates statistically significant trends (p < 0.05) over 3.7%, 7.9%, 77.9%, and 22.5% of the study area, respectively. August therefore remains the month with the strongest and most spatially robust warming signal.
Terra observations were nevertheless retained for the early descriptive part of the 2000–2025 MODIS record. To assess whether the two initial Terra-only years materially influence the spatial characterization of LST extremes, we additionally examined the year in which the absolute maximum LST occurred at each grid cell. Less than 1% of the analyzed area recorded its absolute maximum LST in 2000, while no pixels recorded their absolute maximum in 2001. This indicates that retaining these initial years has only a minor influence on the spatial pattern of extreme LST, and this additional result has been included in the revised manuscript.
We nevertheless acknowledge that the complete descriptive 2000–2025 MODIS record cannot be regarded as a fully homogenized multi-platform climate data record because of the different Terra and Aqua daytime overpass times. This limitation is now explicitly stated in the revised manuscript. However, the principal long-term trend assessment is no longer affected by this issue because it is based exclusively on Aqua observations. The purpose of retaining the early Terra observations is to preserve information from the beginning of the MODIS archive and to provide a complete descriptive assessment of summer surface thermal conditions over the study period, rather than to treat Terra and Aqua observations as fully interchangeable measurements.
We appreciate the reviewer’s reference to Good et al. (2022), which emphasizes the importance of sensor stability and platform consistency in satellite-based climate trend analyses. We agree with this general principle and have consequently adopted a more conservative interpretation of the trends in the revised manuscript. In particular, the revised text reports robust Sen’s slope estimates together with Mann–Kendall significance and explicitly documents the Terra–Aqua platform transition and the different daytime observation times.
Physical distinction between LST, air temperature, and heat stress
Response: We thank the reviewer for highlighting the need to distinguish more clearly between radiometric land surface temperature (LST), near-surface air temperature, and human heat stress. We agree that these quantities represent different physical processes and should not be interpreted interchangeably.
We would first like to clarify that the manuscript does not assume a constant LST–Tmax offset of 2.72 °C. This value represents the mean difference observed in our paired dataset and was reported as a descriptive statistic rather than as a fixed physical relationship between surface and air temperature. We have revised the corresponding text to make this distinction explicit. The LST–Tmax relationship is now presented together with its variability, recognizing that surface–air temperature differences depend on surface characteristics and environmental conditions.
We have also revised the terminology associated with the LST classes. The thresholds of 35, 40, 45, 50, and 55 °C refer exclusively to MODIS radiometric land surface temperature and are used to characterize different levels of surface heating. They are not interpreted as equivalent thresholds of 2-m air temperature or as physiological heat-stress thresholds. Accordingly, references that could imply equivalence between LST, air temperature, and human thermal stress have been removed or reformulated throughout the manuscript.
Finally, we agree that cloud-related sampling must be considered when interpreting MODIS daytime LST. We have therefore removed the previous statement that clear-sky sampling has a “negligible structural impact”. MODIS thermal observations are available primarily under clear-sky conditions, and cloud cover can reduce the spatial and temporal availability of valid LST observations. The revised manuscript explicitly interprets the MODIS results as representing clear-sky surface thermal conditions and acknowledges clear-sky sampling bias as a limitation, particularly during partly cloudy heat-warning days.1. Structure and weighting of the Heat Vulnerability Index (HVI)
Response: We thank the Reviewer for this detailed comment and agree that the conceptual structure and weighting of the HVI required substantial revision and clearer justification. In response to this comment, we redesigned the HVI to more explicitly distinguish the principal dimensions of heat vulnerability and to incorporate adaptive capacity directly into the index.
The revised HVI consists of three equally weighted conceptual dimensions: thermal hazard, demographic sensitivity, and lack of adaptive capacity. Thermal hazard is calculated as the mean of normalized MODIS LST and cumulative heat-warning severity. Demographic sensitivity is represented by the normalized proportion of residents aged ≥65 years within each 1-km grid cell. Adaptive capacity incorporates tree-cover density, accessibility to emergency healthcare, household air-conditioning prevalence, and housing thermal insulation. Because higher adaptive capacity reduces vulnerability, the resulting Adaptive Capacity (AC) component is inverted to obtain Lack of Adaptive Capacity (LAC). The final index is therefore calculated as:
HVI = (Hazard + Sensitivity + LAC) / 3
The weighting structure was also revised to avoid unintended overrepresentation of dimensions for which multiple indicators were available. Within the hazard domain, LST and warning severity receive equal weights. Within adaptive capacity, air-conditioning prevalence and thermal insulation are first combined into a single housing adaptive-capacity component. This housing component is then equally weighted with tree-cover density and healthcare accessibility. Consequently, the three main HVI dimensions—Hazard, Sensitivity, and LAC—each contribute one-third to the final index, while no individual adaptive-capacity indicator receives additional weight simply because more variables were available for a particular domain.
Population exposure is no longer included as a separate component of the HVI. Instead, we introduced a separate Heat Exposure Index (HEI), calculated from normalized LST, cumulative warning severity, and population density. This distinction prevents population concentration from being counted as an additional vulnerability dimension in the HVI and provides a clearer conceptual separation between population exposure and multidimensional vulnerability. Total population is subsequently used only to quantify the number and proportion of residents located within each HVI class and does not contribute to the HVI score itself.
This revision also changes the interpretation of highly populated urban areas such as Bucharest. High population density can contribute to high HEI values, reflecting concentrated population exposure under elevated thermal conditions, but it does not automatically produce the highest HVI values. The revised HVI additionally depends on demographic sensitivity and adaptive capacity. Consequently, the highest vulnerability levels are predominantly concentrated in southern Romania, where elevated thermal hazard coincides with a higher proportion of older residents and comparatively limited adaptive capacity. The HEI and HVI are therefore interpreted as complementary rather than interchangeable indicators of heat exposure and heat vulnerability.
We also addressed the spatial-scale issue associated with the adaptive-capacity indicators. Tree-cover density and healthcare accessibility are represented at the 1-km grid-cell level, whereas household air-conditioning prevalence and thermal-insulation data are available at county level. These county-level values were assigned to the corresponding 1-km grid cells without spatial interpolation and retain their county-level representativeness. This limitation is now explicitly acknowledged, and no artificial within-county spatial variability is inferred from these variables.
Regarding normalization, all HVI components were transformed to a common 0–1 scale using robust min–max normalization based on the 2nd and 98th percentiles. Values below the 2nd percentile were assigned 0 and values above the 98th percentile were assigned 1, thereby limiting the influence of extreme outliers while retaining the relative spatial variation of the indicators. The revised Methods now describe this procedure explicitly.
Finally, we revised the HVI classification to cover the complete theoretical 0–1 range of the index. The previous empirical class limits were replaced by five equal-width classes: Very Low (0.00–<0.20), Low (0.20–<0.40), Moderate (0.40–<0.60), High (0.60–<0.80), and Very High (0.80–1.00). Thus, all possible HVI values are now explicitly covered by the classification scheme. In the revised analysis, HVI values range from approximately 0.081 to 0.929, confirming that the Very High class (>0.80) is represented in the dataset.Circular validation
Response We thank the Reviewer for this important clarification and agree that the relationships reported in the previous version between the HVI and its constituent variables did not constitute independent validation. These relationships are therefore no longer described as validation and are interpreted only as reflecting the internal behavior of the composite index.
In response to the Reviewer’s concern, we additionally introduced an exploratory external health-outcome assessment using independent circulatory-system mortality data from the Romanian National Institute of Statistics (INS). The population-weighted county-level HVI showed a significant positive association with circulatory-system mortality (Spearman ρ = 0.664, p < 0.001, n = 42), and the association remained significant after controlling for the county-level proportion of residents aged ≥65 years (partial Spearman ρ = 0.533, p < 0.001).
Because the available mortality data are annual and county-level rather than heatwave-specific and temporally resolved, this assessment is not interpreted as direct epidemiological validation or evidence of causality. The HVI is therefore presented primarily as a spatial screening and prioritization index, while the mortality analysis provides an exploratory independent assessment of its external health relevance.2. Statistical trend testing and serial correlation
Response: We thank the Reviewer for this important comment. In response, we substantially revised the statistical analysis of the MODIS LST trends. The ordinary least-squares (OLS) analysis used in the previous version of the manuscript was replaced by a non-parametric approach based on Sen’s slope estimator, while the statistical significance of pixel-level monotonic trends was assessed using the Mann–Kendall test. Sen’s slopes are expressed in °C decade⁻¹, and Fig. 5 has been revised to distinguish explicitly between trend magnitude and statistical significance, with hatching identifying pixels with Mann–Kendall p < 0.05.
To eliminate the temporal inhomogeneity associated with the Terra–Aqua platform transition from the trend assessment, the revised trend analysis is based exclusively on Aqua MODIS observations. July and August trends were calculated for 2002–2025 (n = 24 annual observations), whereas June and the complete JJA series were calculated for 2003–2025 (n = 23 annual observations), when a complete Aqua summer record is available. Thus, the trend analysis no longer combines Terra and Aqua observations.
The revised Aqua-only analysis yields mean Sen’s slopes of 0.17 °C decade⁻¹ for June, 0.83 °C decade⁻¹ for July, 2.26 °C decade⁻¹ for August, and 1.10 °C decade⁻¹ for JJA. Positive trends occur over approximately 65%, 88%, 100%, and 98% of the study area, respectively. However, statistically significant Mann–Kendall trends (p < 0.05) occur over only 3.7%, 7.9%, 77.9%, and 22.5% of the respective areas. The revised manuscript therefore clearly distinguishes the spatial prevalence of positive slopes from the considerably smaller spatial extent of statistically significant trends, particularly for June, July, and JJA.
We agree with the Reviewer that the relatively short satellite record and substantial interannual variability require caution in interpreting these trends. Sen’s slope was selected because it provides a robust non-parametric estimate of monotonic trend magnitude and is less sensitive to individual extreme observations than OLS regression. We have therefore revised the interpretation throughout the manuscript to avoid treating all positive slopes as evidence of statistically significant warming.
We also agree that pixelwise significance does not by itself establish field-wide significance. The Mann–Kendall p-values reported in the present analysis represent pixel-level significance and were not interpreted as formal field-wide significance. Finally, the MODIS results are now compared quantitatively with the longer station-based air-temperature record rather than being described as independent “confirmation” of the satellite trends. We emphasize that LST and near-surface air temperature represent different physical quantities and that the station record is therefore used as an independent climatic context for interpreting the satellite-observed surface thermal changes, rather than as a direct validation of the MODIS trend magnitude.
Figure 5v1. Sen’s slope trends (°C decade⁻¹) in summer LST across Romania (2000–2025). Hatched areas indicate Mann–Kendall significance (p < 0.05).
Figure 5v2. Sen’s slope trends (°C decade⁻¹) in summer Aqua MODIS LST across Romania. July and August trends cover 2002–2025 (n = 24), whereas June and JJA trends cover 2003–2025 (n = 23). Hatched areas indicate Mann–Kendall significance (p < 0.05).Finally, the relationship with the independent near-surface air-temperature observations has been reformulated and quantitatively evaluated. Rather than describing the air-temperature record as simply “confirming” the satellite trends, the revised manuscript compares spatially averaged Aqua MODIS daytime LST with daily maximum near-surface air temperature (Tmax) using 2,162 paired country-days from the common observation period (July–August 2002 and JJA 2003–2025). The comparison yielded R² = 0.69, a mean bias of 2.72 °C, and an RMSE of 3.48 °C. These results indicate substantial covariability between the two temperature measures under common large-scale thermal conditions, while also demonstrating systematic differences in their absolute values. We therefore explicitly recognize that satellite-derived LST and near-surface air temperature are physically distinct quantities and that LST should not be interpreted as a direct substitute for air temperature or human heat stress.
6. Missing regional literature
Response: We thank the reviewer for highlighting these important contributions to the Romanian heat-climate and heat-vulnerability literature. We agree that the previous version of the manuscript did not sufficiently acknowledge several relevant national and regional studies, and we have substantially expanded the literature review and Discussion accordingly.
In particular, we have incorporated the country-scale MODIS analysis of surface urban heat-island conditions in Romania by Cheval et al. (2022) and the assessment of heat hazard-risk across 77 Romanian urban areas by Cheval et al. (2023). These studies provide important context for interpreting the spatial patterns of MODIS LST, population exposure, urban structure, and heat-related risk identified in the present analysis. We have also included the earlier MODIS-based analysis of the Bucharest surface urban heat island by Cheval and Dumitrescu (2015), the assessment of heatwave effects on the surface urban heat island in Cluj-Napoca by Herbel et al. (2018), and the analysis of land-use/land-cover change and surface urban heat-island dynamics in Bucharest by Grigoraș and Urițescu (2019).
The revised manuscript also incorporates the regional heat-health vulnerability analysis of Mocanu et al. (2021) and the health-impact studies of Scripca et al. (2022). These studies are particularly relevant for placing our HVI results within the broader Romanian literature on population sensitivity, adaptive capacity, heat-related mortality, and cardiovascular vulnerability. We further acknowledge the recent remote-sensing analysis of heatwave–urban vegetation interactions in Bucharest by Zoran et al. (2026), which provides additional context for MODIS-era heat conditions in the capital.
In response to the reviewer’s comment, we have also revised the statement regarding the novelty of the present study. We do not intend to imply that previous national-scale or regional heat-risk, urban-climate, or heat-vulnerability assessments have not been conducted in Romania. Rather, the specific contribution of the present study lies in the integration of a 26-year MODIS LST record (2000–2025), the long-term record of official operational heat-warning severity, and 1-km population-grid indicators within a spatially continuous framework covering both urban and rural areas across Romania. The revised manuscript now positions this contribution explicitly in relation to the existing Romanian literature rather than presenting it as the first heat-vulnerability analysis of the country.Specific and Minor Comments
L. 153–155: The statement that LST measures the "atmospheric greenhouse effect" is physically incorrect. Please rephrase in terms of radiometric surface temperature and thermal infrared emissions.
Response: Thank you for this observation. We agree that the original reference to the greenhouse effect was physically imprecise. We have revised the sentence to describe LST in terms of surface energy-balance conditions and land–atmosphere interactions.
Figure 1: This overview figure is presented before the Data and Methods section without explaining the underlying data source, spatial aggregation, or statistical regression.
Response: Thank you for this observation. We agree that the original placement of Figure 1 preceded the description of its underlying dataset and analytical procedure. We have therefore moved Figure 1 from the Introduction to the Data and Methods section, after the description of the underlying datasets and analytical workflow. We have also expanded the Data and Methods sections to specify the source and spatial resolution of the ANM Tmax dataset, the June–August temporal aggregation, the spatial averaging across the 13 elevation bands.
Figure 4: Please clearly specify the sample unit for n=2343 (are these station-days, county-days, or pixel-days?).
Response: Thank you for pointing out this ambiguity. We have clarified that the revised analysis includes n = 2,162 paired country-day observations. Each observation consists of the spatially averaged MODIS LST and gridded Tmax over Romania for one summer day with both variables available; it does not represent an individual station, county, or pixel. We have added this information to the Methods section and the Figure caption. We have also corrected the figure title to indicate that the analysis includes summer days with paired observations, rather than exclusively heatwave days.
Figures 3, 5, 6, 10, 11: The text and label sizes across multi-panel figures are much too small to read comfortably. Please enlarge all axis titles, tick labels, and legends.
Response: We thank the reviewer for this comment. Figures 3, 5, 6, 10, and 11 have been carefully reviewed for readability. Text elements, including axis titles, tick labels, panel labels, and legends, have been enlarged where necessary to improve readability and ensure consistent presentation across the multi-panel figures.
L. 376–377: The definition of heatwave days in 2024 (56 days) should be reconciled with the monthly sum of warning days (15 + 20 + 22 = 57 days).
Response: Thank you for identifying this inconsistency. The value of 20 warning days reported for July 2024 was a typographical error. The correct value is 19 days; therefore, the monthly total is 15 + 19 + 22 = 56 heatwave-warning days, consistent with the annual value reported in the manuscript. We have corrected this value in the revised manuscript.
Wording / Style: Please replace informal or promotional phrases such as "empirical blueprint", "actionable diagnostic tool", and "transcends theoretical mapping" with sober scientific descriptions.
Response: Thank you for this observation. We agree that several expressions in the original manuscript were overly promotional and could overstate the demonstrated applicability of the proposed index. We have replaced phrases such as “empirical blueprint”, “actionable diagnostic tool”, and “transcends theoretical mapping” with more neutral scientific descriptions. We have also revised similar wording in the Abstract, Discussion, Future Work, and Conclusions and clarified that the indices provide spatial screening and prioritization tools rather than validated epidemiological predictors of heat-related health outcomes.
Data and Code Availability: In line with Copernicus data policies, please deposit the processing scripts (GEE/Python/R) and processed figure-generating datasets in a permanent FAIR repository (e.g., Zenodo).
Response: We thank the reviewer for this recommendation and agree that transparency and reproducibility are important aspects of the study. We have revised the Data and Code Availability section to provide a clearer description of the accessibility of the datasets and processing materials used in the analysis.
The satellite products used in this study are publicly available from their respective data providers and are identified and cited in the revised manuscript. In contrast, the meteorological warnings provided by the Romanian National Meteorological Administration are subject to institutional data-use and redistribution restrictions and therefore cannot be deposited by the authors in an unrestricted public repository. We also acknowledge the reviewer’s recommendation concerning the processing scripts. The GEE, Python, and R scripts were developed as part of the analytical workflow of the present study. Copernicus Publications encourages the deposition of software, algorithms, and model code in FAIR-aligned repositories whenever possible. At present, the complete processing workflow has not been deposited in a permanent public repository. The revised manuscript therefore provides detailed information on the datasets, processing procedures, statistical methods, and index calculations required to understand and reproduce the analytical approach.
We have consequently revised the Data and Code Availability statement to distinguish between publicly accessible source data, institutionally restricted meteorological data, and author-developed processing code. The latter may be made available by the corresponding author upon reasonable request, subject to applicable institutional requirements.
The Conclusion is repetitive and substantially overstates operational and public-health applicability. Shorten it and include limitations and uncertainty alongside the results.
Response: We thank the reviewer for this comment. The Conclusions section has been substantially revised and shortened to reduce repetition of results and to provide a more balanced interpretation of the study findings. We have also revised statements concerning the operational and public-health applicability of the proposed framework. The Conclusions therefore now emphasize the principal scientific findings while presenting the potential application of the framework more cautiously, as a tool for identifying areas where combined thermal hazard, population exposure, and demographic susceptibility may warrant further investigation and adaptation planning.Minimum requirements for a credible new submission
1. Quantify the potential Terra-Aqua/observation-time discontinuity and incomplete 2002 season using Terra-only, Aqua-only, common-period, exclusion, and breakpoint tests; rebuild or restrict the series if the results show material bias.
Response: We thank the reviewer for emphasizing the Terra–Aqua platform transition and the different daytime observation times in the long-term LST record. We agree that this distinction should be clearly acknowledged when describing and interpreting a multi-decadal satellite temperature record. We would first like to clarify that the summer 2002 record used in the present analysis is not incomplete. Terra MOD11A1 observations were used for June 2002, while Aqua MYD11A1 observations were used from July 2002 onward. This has now been stated explicitly in the revised manuscript.
The objective of the present study, however, is not to construct a fully homogenized satellite climate data record or to quantify platform-specific differences in the manner of a dedicated sensor-intercomparison study. Rather, the MODIS record is used to characterize the spatial and temporal evolution of summer surface thermal conditions across Romania over the 2000–2025 period. For this purpose, the trend methodology has been substantially strengthened in the revised manuscript by replacing ordinary least-squares regression with the non-parametric Sen’s slope estimator, while trend significance is evaluated using the Mann–Kendall test. This approach reduces sensitivity to individual extreme years and provides a more robust estimate of the monotonic trend than the original OLS analysis.
Terra and Aqua have different daytime overpass times, and the transition between the two platforms is now explicitly documented in the revised manuscript. This distinction is considered when interpreting the long-term MODIS record.
Importantly, the principal conclusions of the study do not depend on a single endpoint comparison or on the magnitude of one trend estimate. The revised Aqua-only analysis shows spatially extensive positive Sen’s slopes, particularly in August, when positive trends occur over 100% of the analysed area and statistically significant Mann–Kendall trends (p < 0.05) occur over 77.9%. The JJA mean similarly shows positive trends over approximately 98% of the analysed area, although statistically significant trends are more spatially restricted, occurring over 22.5%. These results therefore distinguish between the widespread occurrence of positive trend estimates and the more limited areas where trends are statistically significant. Together with the quantitative comparison between MODIS daytime LST and near-surface Tmax, these spatial patterns are consistent with a broad intensification of summer surface thermal conditions, while recognizing the relatively short satellite record and the physical differences between LST and near-surface air temperature.
We therefore retained the full 2000–2025 Terra–Aqua record only for descriptive analyses, while restricting the principal trend analysis to the temporally consistent Aqua record (July–August 2002–2025 and June/JJA 2003–2025). This preserves the early MODIS observations without allowing the Terra–Aqua transition to influence the principal trend estimates.2. Fully document QA, cloud/missing-data handling, aggregation, spatial processing, and trend inference.
Response: We thank the reviewer for this recommendation. The Methods section has been revised and expanded to provide a more complete and reproducible description of the MODIS LST processing workflow and statistical analysis.
Specifically, we now provide additional information on the MODIS Collection 6.1 products used, the treatment of quality-control information and invalid observations, the handling of missing and cloud-affected pixels, and the clear-sky nature of the resulting MODIS LST observations. We also clarify the temporal aggregation from daily observations to monthly summer composites and the spatial processing applied to obtain a consistent national-scale raster dataset.
The trend-analysis procedure has also been documented more explicitly. Trends are calculated independently at the pixel level using Aqua MODIS observations for July and August 2002–2025 (n = 24) and for June and JJA 2003–2025 (n = 23). The original ordinary least-squares approach has been replaced by Sen’s slope estimator, expressed in °C decade⁻¹, while statistical significance is evaluated using the Mann–Kendall test at p < 0.05. Figure 5 has been revised accordingly to distinguish trend magnitude and direction from statistical significance.
Finally, the revised manuscript explicitly acknowledges that MODIS thermal-infrared LST observations are primarily available under clear-sky conditions. Cloud-affected or otherwise invalid observations are treated as missing rather than being spatially interpolated, and no attempt is made to reconstruct LST beneath clouds. This limitation is now stated explicitly when interpreting the resulting long-term surface-temperature patterns.3. Separate LST exceedances, air-temperature heatwaves, and operational warnings.
Response: We thank the reviewer for this important recommendation. We agree that land-surface temperature exceedances, air-temperature-based heatwave conditions, and operational heat warnings represent distinct quantities and should not be used interchangeably.
The revised manuscript therefore distinguishes these three components explicitly throughout the Methods, Results, figure captions, and Discussion. MODIS LST exceedances refer exclusively to radiometric land-surface temperature and are used to characterize the magnitude, spatial extent, and persistence of surface heating. The LST thresholds used in the analysis are consequently described as surface-temperature thresholds and are not interpreted as equivalent thresholds of near-surface air temperature or human heat stress.
Near-surface air temperature (Tmax) is treated separately as an atmospheric temperature variable. Its relationship with MODIS daytime LST was evaluated independently using 2,162 paired country-days from the common Aqua MODIS–Tmax observation period (July–August 2002 and JJA 2003–2025). The comparison yielded R² = 0.69, a mean bias of 2.72 °C, and an RMSE of 3.48 °C. This analysis was used to quantify the relationship between surface and near-surface air temperature rather than to convert LST into air temperature or to assume a fixed LST–Tmax offset. LST and Tmax were therefore treated as physically distinct temperature variables throughout the analysis.
Operational heat warnings issued by the Romanian National Meteorological Administration are likewise treated as a separate information source. They represent operational assessments of hazardous heat conditions and are used to characterize the occurrence, duration, and severity of officially identified heat episodes. They are not interpreted as direct equivalents of either MODIS LST thresholds or a single air-temperature threshold.
We have revised the terminology throughout the manuscript accordingly and now explicitly distinguish between surface thermal conditions identified from MODIS LST, near-surface air-temperature conditions, and operational heat-warning information. This separation is also reflected in the interpretation of the HVI, in which LST and operational warning severity represent distinct components of the thermal-hazard dimension.4. Restrict warning-trend conclusions to a homogeneous period or provide validated harmonization.
Response: We thank the reviewer for this important comment. The long-term analysis combines the available operational heat-warning information with an air-temperature-based criterion for the earlier period, when comparable operational warning records were not available. For this earlier period, a daily maximum air-temperature threshold of 35°C was used to identify severe hot conditions, while the subsequent period was characterized using the official heat-warning information issued by the Romanian National Meteorological Administration. This approach was adopted to provide temporal continuity in the characterization of major summer heat episodes across the 2000–2025 study period. We have clarified the distinction between the two sources of information in the manuscript and interpret the resulting long-term record accordingly.3. Redesign and rename the HVI/risk index; justify indicators, weights, normalization, and classes; quantify sensitivity and uncertainty.
Response: We thank the Reviewer for this important comment. We agree that the conceptual structure, terminology, indicator selection, weighting scheme, normalization procedure, class definition, and sensitivity of the index required clearer justification. In response, we substantially redesigned the vulnerability framework and revised the corresponding Methods, Results, Discussion, and Conclusions.
The revised Heat Vulnerability Index (HVI) is structured around three complementary and equally weighted conceptual domains: thermal hazard, demographic sensitivity, and lack of adaptive capacity. Thermal hazard combines normalized MODIS LST and cumulative heat-warning severity; demographic sensitivity is represented by the normalized proportion of residents aged ≥65 years; and adaptive capacity incorporates tree-cover density, accessibility to emergency healthcare, household air-conditioning prevalence, and housing thermal insulation. Because greater adaptive capacity reduces vulnerability, the Adaptive Capacity (AC) component is inverted to obtain Lack of Adaptive Capacity (LAC). The final index is therefore calculated as:
HVI = (Hazard + Sensitivity + LAC) / 3
Following this redesign, we retained the term Heat Vulnerability Index (HVI), because the revised formulation explicitly represents thermal hazard, demographic sensitivity, and adaptive capacity, while population exposure is treated separately rather than being incorporated as an additional vulnerability component.
The weighting structure was defined at the conceptual-domain level rather than by assigning equal weights to all individual indicators. Within the hazard domain, LST and heat-warning severity are equally weighted. Within adaptive capacity, air-conditioning prevalence and thermal insulation are first combined into a single housing adaptive-capacity component; this housing component is then equally weighted with tree-cover density and healthcare accessibility. This domain-balanced structure prevents the housing dimension from receiving disproportionate influence simply because it is represented by two indicators. Hazard, Sensitivity, and LAC consequently each contribute one-third to the final HVI. Equal domain weights were selected as a transparent baseline in the absence of sufficient empirical evidence to justify preferential weighting of one vulnerability dimension over another.
Population exposure was separated from the HVI to avoid double representation of population characteristics. A complementary Heat Exposure Index (HEI) was introduced to characterize the spatial coincidence of normalized LST, cumulative heat-warning severity, and normalized population density. Total population does not enter the HVI calculation and is used subsequently only to quantify the number and proportion of residents located within each vulnerability class. This revision provides a clearer distinction between heat exposure and multidimensional heat vulnerability.
To further preserve the distinction between exposure, demographic sensitivity, and adaptive capacity, we introduced a complementary Heat Vulnerability Prioritization (HVP) framework. Rather than combining these dimensions into a single continuous score, HVP independently classifies HEI, demographic sensitivity, and adaptive capacity into distribution-based tertiles, yielding 27 possible profiles. The High HEI–High Sensitivity–Low Adaptive Capacity (H–H–L) profile identifies areas where all three adverse conditions coincide. For communication and spatial screening, the 27 profiles were additionally summarized into Highest Priority, Elevated Priority, and Other combinations according to the number of adverse conditions present. The Highest Priority class comprises 3,144 1-km grid cells and 593,581 residents, corresponding to 3.11% of the national population. HVP therefore complements, rather than replaces, the continuous HVI by explicitly identifying locations where adverse exposure, demographic sensitivity, and adaptive-capacity conditions coincide.
All component indicators were transformed to a common 0–1 scale using robust min–max normalization based on the 2nd and 98th percentiles. Values below the 2nd percentile were clipped to 0 and values above the 98th percentile to 1. This procedure reduces the influence of extreme observations while retaining spatial variation across the large majority of the distribution. The normalization procedure and clipping rules are now explicitly described in the revised Methods.
The HVI classification was also redesigned to cover the complete theoretical 0–1 range using five equal-width classes: Very Low (0.00–<0.20), Low (0.20–<0.40), Moderate (0.40–<0.60), High (0.60–<0.80), and Very High (0.80–1.00). This replaces the previous empirically restricted classification and ensures that every possible HVI value is assigned to a class. The revised HVI ranges from approximately 0.081 to 0.929, and the Very High class is therefore represented in the final dataset.
To quantify sensitivity to the weighting assumptions, we additionally performed a weighting sensitivity analysis. The baseline equal-domain weighting (Hazard = 1/3, Sensitivity = 1/3, LAC = 1/3) was compared with three alternative scenarios emphasizing Hazard (0.40/0.30/0.30), Sensitivity (0.30/0.40/0.30), or LAC (0.30/0.30/0.40). The resulting spatial rankings remained highly consistent with the baseline HVI, with Spearman rank correlations of 0.9916, 0.9911, and 0.9938, respectively. The spatial overlap of High and Very High vulnerability cells, quantified using the Jaccard index, was 92.02%, 88.33%, and 91.92%, respectively, while county-level rankings were also highly stable, with Spearman correlations of approximately 0.99. These results indicate that the principal spatial patterns are robust to moderate changes in the relative weighting of the three HVI domains.
We nevertheless recognize that this sensitivity analysis evaluates robustness to the tested weighting assumptions rather than demonstrating that the selected formulation is uniquely optimal. Alternative indicator selections, normalization methods, or substantially different weighting assumptions could produce different spatial patterns. Spatial-scale uncertainty also remains because tree-cover density and healthcare accessibility are represented at grid-cell level, whereas air-conditioning prevalence and thermal-insulation indicators are available at county level. These limitations are now explicitly acknowledged in the revised manuscript.
Finally, an exploratory external health-outcome assessment was added using independent county-level circulatory-system mortality data from the Romanian National Institute of Statistics (INS). The population-weighted county-level HVI showed a significant positive association with mean annual circulatory-system mortality (Spearman ρ = 0.664, p < 0.001, n = 42), and the association remained significant after controlling for the county-level proportion of residents aged ≥65 years (partial Spearman ρ = 0.533, p < 0.001). Because the available mortality data are annual and county-level rather than temporally resolved heatwave-specific health records, this analysis is interpreted as an exploratory external health-outcome assessment rather than direct epidemiological validation. Accordingly, the revised manuscript presents the HVI as a transparent and reproducible spatial screening and prioritization index rather than as a validated predictor of heat-attributable mortality or individual-level risk. More rigorous evaluation using temporally resolved, cause-specific health-impact data remains an important direction for future research4. Validate against independent health or impact outcomes, or explicitly present the output as an unvalidated screening index.
Response: We thank the Reviewer for this important comment and agree that evaluation against independent health-impact data is important for assessing the external relevance of the proposed vulnerability framework. In response to this comment, we added an exploratory external health-outcome assessment using independent circulatory-system mortality data from the Romanian National Institute of Statistics (INS).
Because the HVI was developed at 1-km resolution whereas the available mortality data were reported at county level, the 1-km HVI values were aggregated to county level using population-weighted means and compared with mean annual circulatory-system mortality rates derived from the available annual INS data for 2000–2025. The association was evaluated using Spearman rank correlation. Because the proportion of residents aged ≥65 years is included in the HVI as the demographic-sensitivity component and is also intrinsically related to mortality, we additionally performed a partial Spearman correlation controlling for the county-level elderly population share.
The population-weighted county-level HVI showed a significant positive association with circulatory-system mortality (Spearman ρ = 0.664, p < 0.001, n = 42). Importantly, the association remained significant after controlling for the proportion of residents aged ≥65 years (partial Spearman ρ = 0.533, p < 0.001), indicating that the observed relationship was not solely attributable to the demographic-sensitivity component of the HVI. The relationship is presented in the revised manuscript and in Fig. S11.
We nevertheless interpret this analysis cautiously. The available mortality data consist of annual county-level totals rather than daily, heatwave-specific, or age-specific mortality records. Consequently, the observed association cannot be interpreted as direct validation of heat-attributable mortality or as evidence of a causal relationship. We therefore continue to present the HVI primarily as a spatial screening and prioritization index rather than as a validated epidemiological predictor. Relationships between the HVI and its constituent components, as well as the alternative-weight sensitivity analyses, are interpreted only as measures of internal behavior and robustness.
The revised Discussion and Conclusions now explicitly acknowledge these limitations and emphasize that future evaluation using temporally resolved, cause-specific health-impact data would provide a more rigorous epidemiological assessment of the proposed indices.5. Archive code and figure-reproduction data; provide reviewer access to restricted inputs.
Response: We thank the reviewer for this recommendation and recognize the importance of transparency and reproducibility. The materials used in this study are subject to different access and redistribution conditions, which are now more clearly distinguished in the revised Data and Code Availability statement.
The original MODIS LST products used in the study are publicly available from their respective data providers. The Google Earth Engine (GEE) code used for the retrieval and processing of the MODIS LST data is being prepared for public archiving and may be deposited in a permanent repository such as Zenodo following final documentation and quality control of the workflow. The corresponding repository information and DOI can subsequently be provided with the final version of the manuscript, where applicable.
In contrast, the operational heat-warning datasets provided by the Romanian National Meteorological Administration are subject to institutional data-use and redistribution restrictions. The authors are therefore not authorized to deposit these datasets in an unrestricted public repository. However, the daily gridded maximum air temperature data used in the study are already publicly available through the Romanian Open Data Portal (https://data.gov.ro/dataset/date-meteorologice-zilnice-gridate).
Regarding figure-reproduction data, openly redistributable derived products associated with the MODIS LST analysis may also be included with the archived reproducibility materials where appropriate. However, datasets containing or derived directly from restricted institutional information cannot be redistributed where this would conflict with the conditions governing the original data.
We have revised the Data and Code Availability statement accordingly to distinguish clearly between publicly available satellite data, reproducibility materials that can be archived, and institutionally restricted meteorological information.6. Reconstruct the literature review and verify every bibliographic record.
Response: We thank the Reviewer for this recommendation. The literature review has been substantially revised and expanded to provide a more complete and balanced context for the study, with particular attention to previous research on MODIS-derived LST, heat exposure, heat vulnerability, adaptive capacity, and heat-related health impacts in Romania and the wider regional context.
In particular, the revised manuscript now incorporates and discusses the relevant Romanian and regional studies highlighted by the Reviewer, including previous work on MODIS-based surface urban heat islands, urban heat hazard and risk, heat-health vulnerability, land-use and vegetation controls on surface temperature, and temperature-related mortality. These studies are now more explicitly related to the objectives, methodological choices, and contribution of the present analysis.
The novelty statements have also been revised to avoid implying that satellite-based heat conditions or heat vulnerability have not previously been investigated in Romania. Instead, the contribution of the present study is positioned more specifically in terms of the long-term, national-scale integration of MODIS LST observations, heat-warning information, 1-km gridded demographic data, and multiple indicators of adaptive capacity within a spatially continuous framework that distinguishes population exposure (HEI) from multidimensional heat vulnerability (HVI).
In addition, the reference list and corresponding in-text citations were systematically reviewed for bibliographic accuracy and consistency. Author names, publication years, article titles, journal names, volume and issue information, page ranges or article numbers, and DOI information were checked and corrected where necessary. In-text citations were also cross-checked against the reference list to identify missing, duplicated, or inconsistent records.
-
AC2: 'Reply on RC1', Anisoara Irimescu, 16 Sep 2026
-
RC2: 'Comment on egusphere-2026-4415', Anonymous Referee #2, 05 Sep 2026
Authors present a national-scale assessment of summer heatwaves and heat vulnerability in Romania using a 26-year archive (2000–2025) of MODIS Land Surface Temperature (LST) integrated with official meteorological warnings and high-resolution population-grid data.
It’s the first study of this kind for Romania, and very timely given the increased frequency and intensity of heat waves in the Country and, at larger scale, in Europe.
Authors used daily satellite observations to analyze long-term trends, anomalies, and consecutive hot days from two different satellites and develop a 1-km Heat Vulnerability Index (HVI) based on hazard, exposure and vulnerability.
Regarding the hazard, Authors use Land Surface Temperature as representative of the environmental conditions, then include, which might represent a welcome an interesting novelty, an index indicating the severity of the event, derived from the warnings given to the population. Not all the warnings have the same weight, but it increases, with the class (yellow, orange, red). None is said about how these warnings are issued. Are they related to the Air Temperature at 2m or to any other variable? Is it necessary to use this information or is it redundant, and the hot spell is already well represented by LST?
- I ask the Authors to clarify this aspect more in detail.
The idea of using two satellites covering two consecutive period Terra Modis 2000-2001 and Aqua Modis 2002-2025, can be, in principle, a good one, but the overpass over Romania is not at the same time of the day, and to add 2 years of data, does not add any to the study, and could increase a bias or produce a misleading result (during hot periods, LST can change of few degrees in 2 hours, so not representing comparable conditions). I suggest:
- To make the analysis only using the AQUA Modis (or AQUA Terra) for the time available.
- Or make two different studies each using a full set offered by each satellite and then compare the results for the overlap period (if I’m correct one of the two started a bit later)
In the introduction I don’t see the purpose to revise all the different definitions of heat waves which are based on air temperature (max over absolute or percentile threshold) since they are not used in the paper. A single citation is enough.
There is a bit of confusion throughout the manuscript between LST and Air Temperature. If Authors use only LST to calculate the HVI, but they intend to offer a more complete overview of the heatwaves characteristics in different regions of the Country, I suggest to make the text more clear for the reader.
HVI calculations: in the formula adopted Authors consider only the normalized proportion of population aged ≥65 years as vulnerability index, which good, but not enough. Following Bao et al (2015) and Quian and Liu (2025) Vulnerability is a function of the character, magnitude, and rate of climate change as well as the variation to which a system is exposed, its sensitivity, and its adaptive capacity. None of the mentioned elements is taken into account in the present study, so result can be somehow interesting, but far for being of importance in identifying populations and areas most susceptible to the adverse impacts of extreme heat.
Few more scientific articles can be of interest and added to the references
Mărculeț, Cătălina, and Cristina Dumitrică. "The excessive heatings in the Romanian Plain." Central European Journal of Geography and Sustainable Development 2.1 (2020): 30-37.
Onțel, Irina, et al. "Influence of environmental factors on land surface temperature and surface urban heat island. A cross-country analysis in Romania." Sustainable Cities and Society 128 (2025): 106454.
Papathoma-Koehle, Maria, et al. "A common methodology for risk assessment and mapping for south-east Europe: an application for heat wave risk in Romania." Natural Hazards 82.Suppl 1 (2016): 89-109.
I suggest to reject the manuscript in its present form, revise it thoroughly following all the comments and submit a fully revised version.
Citation: https://doi.org/10.5194/egusphere-2026-4415-RC2 -
AC1: 'Reply on RC2', Anisoara Irimescu, 16 Sep 2026
We sincerely thank the Reviewer for the careful evaluation of our manuscript and for the constructive comments and suggestions. We have carefully considered all the points raised and revised the manuscript accordingly.
1. Regarding the hazard, Authors use Land Surface Temperature as representative of the environmental conditions, then include, which might represent a welcome an interesting novelty, an index indicating the severity of the event, derived from the warnings given to the population. Not all the warnings have the same weight, but it increases, with the class (yellow, orange, red). None is said about how these warnings are issued. Are they related to the Air Temperature at 2m or to any other variable? Is it necessary to use this information or is it redundant, and the hot spell is already well represented by LST?
Response: We thank the Reviewer for this important comment. This aspect has been clarified in the revised manuscript. The Romanian colour-coded heat-warning system follows the MeteoAlarm/EMMA framework and the national warning methodology. Yellow, orange, and red heat warnings are based on the maximum air temperature at 2 m, the Temperature–Humidity Index (THI), and the persistence of heatwave conditions, with increasing warning levels reflecting increasing meteorological severity.
We consider the warning-severity information complementary rather than redundant with MODIS LST. LST characterizes the radiometric thermal state of the land surface at the satellite overpass time, whereas the operational warnings incorporate atmospheric conditions relevant to human heat exposure, including near-surface air temperature, humidity, and persistence of extreme heat. Thus, LST provides spatial information on surface heating, while the warning-severity component provides information on the meteorological intensity and persistence of the event as evaluated by the national warning system. Their combination was therefore intended to represent complementary dimensions of heat hazard rather than duplicate the same thermal information.
2. The idea of using two satellites covering two consecutive period Terra Modis 2000-2001 and Aqua Modis 2002-2025, can be, in principle, a good one, but the overpass over Romania is not at the same time of the day, and to add 2 years of data, does not add any to the study, and could increase a bias or produce a misleading result (during hot periods, LST can change of few degrees in 2 hours, so not representing comparable conditions).
Response: We thank the Reviewer for raising this important methodological point. We agree that the different daytime overpass times of Terra and Aqua may influence absolute LST values because of the strong diurnal variability of land surface temperature. We would like to clarify that the two platforms were not combined simultaneously: Terra MODIS was used for 2000–2001 and June 2002, whereas Aqua MODIS was used from July 2002 through 2025. The purpose of including the earlier Terra observations was to extend the MODIS record to the beginning of the satellite archive and retain information from the early 2000s.
Following the Reviewer’s recommendation, we revised the trend analysis to avoid the potential bias associated with combining observations acquired at different daytime overpass times. The Sen’s slope and Mann–Kendall analyses presented in the revised manuscript are therefore based exclusively on Aqua MODIS observations (July–August 2002 and JJA 2003–2025). The revised Aqua-based analysis shows mean Sen’s slopes of 0.17, 0.83, and 2.26 °C decade⁻¹ for June, July, and August, respectively, with a JJA mean of 1.10 °C decade⁻¹. Positive trends occur over 65%, 88%, and 100% of Romania in June, July, and August, respectively, and over 98% for JJA; statistically significant trends (Mann–Kendall, p < 0.05) cover 3.7%, 7.9%, 77.9%, and 22.5% of the study area, respectively. Thus, the principal long-term trend assessment is now derived from a temporally consistent Aqua record rather than from the combined Terra–Aqua series.
We nevertheless retained the Terra observations for the early part of the descriptive 2000–2025 record. To assess whether these initial Terra-only years substantially influence the spatial characterization of LST extremes, we additionally examined the year in which the absolute maximum LST occurred at each grid cell. Less than 1% of the analyzed area recorded its 2000–2025 absolute maximum LST in 2000, while no pixels recorded their absolute maximum in 2001. This indicates that retaining the two initial Terra-only years has only a minor influence on the spatial pattern of extreme LST. This additional result has been included in the Supplementary analysis and the Terra–Aqua difference in overpass time is now explicitly acknowledged in the revised manuscript.
Figure 5v1. Sen’s slope trends (°C decade⁻¹) in summer LST across Romania (2000–2025). Hatched areas indicate Mann–Kendall significance (p < 0.05).
Figure 5v2. Sen’s slope trends (°C decade⁻¹) in summer Aqua MODIS LST across Romania. July and August trends cover 2002–2025 (n = 24), whereas June and JJA trends cover 2003–2025 (n = 23). Hatched areas indicate Mann–Kendall significance (p < 0.05).3. In the introduction I don’t see the purpose to revise all the different definitions of heat waves which are based on air temperature (max over absolute or percentile threshold) since they are not used in the paper. A single citation is enough.
Response: We thank the Reviewer for this comment. We agree that the air-temperature-based definitions presented in the Introduction are not directly applied as thresholds to the MODIS LST analysis. However, we have retained this contextual overview because its purpose is to provide the general climatological and operational context of heatwave identification rather than to define the satellite-derived LST thresholds used in the study. In particular, the examples illustrate that there is no universally applicable definition of a heatwave and that the temperature thresholds, duration criteria, and consideration of daytime and nighttime conditions vary according to regional climatic conditions and operational practices. This context is particularly relevant for Romania, where temperature thresholds that may characterize heatwave conditions elsewhere can represent normal summer conditions.
To avoid any possible confusion, we have clarified that these conventional heatwave definitions refer to near-surface air temperature and are presented only as background context, whereas the satellite component of our study specifically analyzes land surface temperature (LST). The two quantities are not treated as interchangeable in the manuscript.
4. There is a bit of confusion throughout the manuscript between LST and Air Temperature. If Authors use only LST to calculate the HVI, but they intend to offer a more complete overview of the heatwaves characteristics in different regions of the Country, I suggest to make the text more clear for the reader.
Response: We thank the Reviewer for this important comment. We agree that a clear distinction between land surface temperature (LST) and near-surface air temperature (Tmax) is essential for the interpretation of the study. We have therefore revised the manuscript throughout to distinguish consistently between surface thermal conditions derived from MODIS LST and atmospheric thermal conditions represented by Tmax.
MODIS LST is used to characterize the spatial and temporal patterns of land-surface heating, whereas Tmax is treated separately as a near-surface atmospheric temperature variable. The relationship between the two variables was quantitatively evaluated using 2,162 paired country-day observations from the common Aqua MODIS–Tmax period (July–August 2002 and JJA 2003–2025), yielding R² = 0.69, a mean bias of 2.72 °C, and an RMSE of 3.48 °C. These results are interpreted as evidence of substantial covariability under common large-scale thermal conditions, while recognizing that LST and Tmax are physically distinct quantities. Accordingly, LST is not interpreted as a direct substitute for near-surface air temperature or human heat stress.
We also clarify that the revised Heat Vulnerability Index (HVI) is not calculated from LST alone. The HVI consists of three equally weighted conceptual domains: thermal hazard, demographic sensitivity, and lack of adaptive capacity. The thermal-hazard component combines normalized MODIS LST and cumulative heat-warning severity; demographic sensitivity is represented by the normalized proportion of residents aged ≥65 years; and lack of adaptive capacity is derived from tree-cover density, accessibility to emergency healthcare, household air-conditioning prevalence, and housing thermal insulation. Population exposure is treated separately through the Heat Exposure Index (HEI), rather than being incorporated into the HVI.
We have consequently revised the terminology in the Methods, Results, figure captions, Discussion, and Conclusions to avoid using LST, air temperature, and heat stress interchangeably. The broader characterization of heatwave conditions in Romania therefore draws on complementary but physically distinct information from MODIS LST, near-surface Tmax, and operational heat-warning records, while the HVI uses LST specifically as one component of its thermal-hazard domain.5. HVI calculations: in the formula adopted Authors consider only the normalized proportion of population aged ≥65 years as vulnerability index, which good, but not enough. Following Bao et al (2015) and Quian and Liu (2025) Vulnerability is a function of the character, magnitude, and rate of climate change as well as the variation to which a system is exposed, its sensitivity, and its adaptive capacity. None of the mentioned elements is taken into account in the present study, so result can be somehow interesting, but far for being of importance in identifying populations and areas most susceptible to the adverse impacts of extreme heat.
Response: We thank the Reviewer for this important comment. We agree that heat vulnerability is multidimensional and should not be represented solely by demographic sensitivity. Following the Reviewer’s recommendation and the conceptual framework highlighted by Bao et al. (2015) and Qian and Liu (2025), we substantially revised the vulnerability framework to explicitly incorporate thermal hazard, demographic sensitivity, and adaptive capacity.
In the revised manuscript, the proportion of residents aged ≥65 years represents only the demographic sensitivity component and no longer constitutes the principal vulnerability information by itself. Thermal hazard is represented by normalized MODIS LST and cumulative heat-warning severity, while adaptive capacity is now explicitly incorporated using four indicators describing environmental, healthcare, and residential capacity to cope with extreme heat: tree-cover density, accessibility to emergency healthcare facilities, household air-conditioning prevalence, and housing thermal insulation.
To avoid disproportionately weighting the residential dimension because two housing indicators were available, air-conditioning prevalence and thermal insulation were first combined into a single housing adaptive-capacity domain. The resulting adaptive-capacity component therefore gives equal weight to tree-cover density, healthcare accessibility, and housing adaptive capacity. Because greater adaptive capacity reduces vulnerability, this component was inverted to represent Lack of Adaptive Capacity (LAC). The revised HVI is consequently calculated as the equal-weighted combination of thermal Hazard, demographic Sensitivity, and Lack of Adaptive Capacity:
HVI = (Hazard + Sensitivity + LAC) / 3
This revision substantially expands the original formulation and directly addresses the Reviewer’s concern that adaptive capacity was not represented in the index. The revised HVI therefore distinguishes three complementary dimensions of vulnerability: the intensity and cumulative occurrence of heat conditions, the demographic sensitivity of the exposed population, and the capacity of local environmental, healthcare, and residential systems to mitigate or respond to heat.
We also recognize that some adaptive-capacity information is not available nationally at 1-km resolution. Tree-cover density and healthcare accessibility are represented at grid-cell level, whereas air-conditioning prevalence and thermal-insulation data are available at county level. Rather than artificially interpolating these variables, county values were assigned to the corresponding 1-km grid cells and explicitly retained as county-level contextual characteristics. The revised manuscript now clearly identifies this scale mismatch as a limitation and does not interpret these indicators as representing within-county household-level variability.
In addition, we introduced a separate Heat Exposure Index (HEI), incorporating LST, cumulative warning severity, and population density, to distinguish population exposure from multidimensional vulnerability. We also developed a complementary Heat Vulnerability Prioritization (HVP) framework that independently combines HEI, demographic sensitivity, and adaptive capacity to identify locations where high exposure, high sensitivity, and low adaptive capacity coincide. These revisions provide a clearer conceptual separation between exposure, sensitivity, and adaptive capacity and considerably strengthen the original vulnerability assessment.6. Few more scientific articles can be of interest and added to the references
Response: We thank the Reviewer for these valuable suggestions. All three recommended studies have been included and cited in the revised manuscript, as they provide relevant complementary information for the present analysis.
We note that Onțel et al. (2025) used Sentinel-3 data for LST analysis in Romania. This is particularly relevant in the context of the Reviewer’s previous observation concerning satellite acquisition times, since the nominal daytime equatorial crossing time of Sentinel-3 SLSTR is approximately 10:00 local solar time, close to that of Terra MODIS (~10:30), whereas Aqua MODIS has a nominal daytime equatorial crossing time of ~13:30 local solar time. This further illustrates that LST observations acquired at different daytime overpass times are used in satellite-based thermal studies, while the acquisition time must be explicitly considered when interpreting and comparing the resulting LST values.
We also remark that retaining the beginning of the MODIS record provides relevant climatological information. The year 2000 represents an important early extreme-heat year in the record: Mărculeț and Dumitrica (2020), one of the studies recommended by the Reviewer, identified 2000 as one of the major record-breaking extreme-heating years in the Romanian Plain, with maximum air temperatures exceeding 42 °C at several meteorological stations. The same study also used MODIS observations to characterize the spatial distribution of LST during the extreme July 2000 episode. Therefore, retaining the early Terra period allows this important extreme event to be represented in our long-term satellite record.
We thank the Reviewer for this important comment and for recommending Papathoma-Koehle et al. (2016). We agree that heat vulnerability is multidimensional and that a comprehensive assessment should ideally consider exposure, sensitivity, and adaptive capacity.
We would like to clarify that, in the revised framework, the proportion of residents aged ≥65 years represents specifically the demographic sensitivity component rather than vulnerability as a whole. In response to the Reviewer’s comment, the HVI formulation was substantially revised to explicitly incorporate three complementary dimensions: thermal hazard, demographic sensitivity, and adaptive capacity. Thermal hazard combines normalized MODIS LST and cumulative heat-warning severity, demographic sensitivity is represented by the normalized proportion of residents aged ≥65 years, and adaptive capacity incorporates tree-cover density, accessibility to emergency healthcare, household air-conditioning prevalence, and housing thermal insulation. Because greater adaptive capacity reduces vulnerability, this component is inverted and expressed as Lack of Adaptive Capacity (LAC). The revised HVI therefore integrates these three dimensions with equal weights:
HVI = (Hazard + Sensitivity + LAC) / 3
Adaptive capacity was explicitly incorporated into the revised HVI using four indicators representing environmental, healthcare, and housing-related dimensions: tree-cover density, accessibility to emergency healthcare, household air-conditioning prevalence, and housing thermal insulation. Tree-cover density and healthcare accessibility were represented at the 1-km grid-cell level, whereas air-conditioning prevalence and thermal-insulation data were available from the 2021 Population and Housing Census at county level. To avoid introducing artificial within-county spatial variability, these county-level values were assigned uniformly to the corresponding 1-km grid cells without spatial interpolation and were therefore interpreted as contextual county-level adaptive-capacity characteristics.
Air-conditioning prevalence and thermal insulation were first combined into a single housing adaptive-capacity domain, which was subsequently integrated with tree-cover density and healthcare accessibility. This domain-balanced structure prevents the housing dimension from receiving disproportionate weight simply because it is represented by two indicators. The resulting Adaptive Capacity (AC) was then inverted to obtain Lack of Adaptive Capacity (LAC), which was incorporated directly into the revised HVI.The study recommended by the Reviewer, Papathoma-Koehle et al. (2016), also emphasizes that detailed heat-vulnerability assessment requires extensive socioeconomic and demographic information that is frequently unavailable, restricted for privacy reasons, or unavailable in an appropriate spatial format. Their Romanian case study therefore adopted alternative indicators where detailed vulnerability data were unavailable. This methodological issue is particularly relevant to our national-scale analysis, for which maintaining spatial consistency at 1-km resolution was a primary consideration.
We have cited Papathoma-Koehle et al. (2016) in the revised Discussion and have clarified this limitation of the present HVI. Accordingly, the proposed HVI is interpreted as a nationally consistent spatial screening index rather than as an exhaustive representation of all dimensions of heat vulnerability.
Data sets
MODIS/Terra MOD11A1 Version 6.1 Z. Wan et al. https://doi.org/10.5067/MODIS/MOD11A1.061
MODIS/Aqua MYD11A1 Version 6.1 Z. Wan et al. https://doi.org/10.5067/MODIS/MYD11A1.061
Eurostat Census 2021 population grid M. Tucci et al. https://ec.europa.eu/eurostat/web/gisco/geodata/population-distribution/population-grids
NMA daily gridded maximum temperature data NMA https://data.gov.ro/dataset/date-meteorologice-zilnice-gridate
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 125 | 57 | 20 | 202 | 36 | 23 | 15 |
- HTML: 125
- PDF: 57
- XML: 20
- Total: 202
- Supplement: 36
- BibTeX: 23
- EndNote: 15
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
General Impression and Summary
The authors present an analysis combining MODIS Land Surface Temperature (LST) data, operational meteorological warning archives from the Romanian National Meteorological Administration (ANM), and 1-km gridded census data to evaluate heatwave trends and spatial heat vulnerability across Romania.
The topic is certainly very timely and relevant for NHESS. Southeastern Europe and the Lower Danube basin are well-known climate change hotspots, and high-resolution spatial assessments combining hazard and demographic data are needed for regional adaptation planning. The attempt to incorporate official meteorological warning archives alongside satellite observations is in principle interesting.
However, after a thorough reading of the text and methods, I have substantial methodological and conceptual concerns that undermine the reliability of the main findings. Most critically, the 26-year trend analysis splices Terra and Aqua MODIS observations without accounting for the ~3-hour difference in overpass times, creating an artificial warming artifact that appears to explain most of the calculated trend. Furthermore, the Heat Vulnerability Index (HVI) is structurally dominated by population density (a 2:1 weighting over hazard), lacks adaptive capacity, and is validated in a circular manner against its own input variables. In addition, there are several problematic entries in the reference list (including non-resolving DOIs and unrelated citations) and several key Romanian urban climate studies have been overlooked.
Given the extent of the recalculations and structural revisions required, I cannot recommend the manuscript for publication in its current form. Below I detail my major concerns and several specific points that the authors should address in a prospective new submission.
Major Concerns
1. Inhomogeneous satellite time series (Terra vs. Aqua overpass times) and 2002 data gap
The most severe methodological issue lies in the construction of the 2000–2025 LST time series (Section 2.1, L. 164–180). The authors state that they used Terra MOD11A1 for 2000–2001 and Aqua MYD11A1 for 2002–2025.
Terra and Aqua have nominal daytime overpass times of approximately 10:30 and 13:30 local solar time, respectively. Consequently, appending Terra observations for 2000–2001 to an Aqua series for 2002–2025 may introduce an observation-time discontinuity. Under the simplifying assumption of a constant Aqua–Terra daytime LST offset of 2.5–4.0 K, the induced ordinary least-squares trend over 2000–2025 would be approximately +0.41 to +0.66 °C/decade. This is potentially substantial relative to the reported JJA trend of +0.92 °C/decade, but the actual bias cannot be established through an assumed offset. It should be quantified using spatially and temporally matched Terra–Aqua observations during their common operating period, stratified by month, land cover, and elevation. Terra-only, Aqua-only, common-period, and platform-transition sensitivity analyses are required before the reported trend magnitude can be considered reliable.
Furthermore, Aqua MYD11A1 (Collection 6.1) data only begin on 4 July 2002. Consequently, June 2002 and early July (over 30 days) are completely missing from the Aqua record. It is unclear how the authors computed full JJA summer statistics for 2002 under these circumstances.
To make the trend analysis defensible, the authors must:
The fact that MODIS has been used successfully for climate-trend studies does not remove this concern. For example, Good et al. (2022) (https://doi.org/10.1029/2022EA002317) evaluated MODIS Terra and Aqua trends separately over the common Aqua period, used monthly anomalies and station-collocated observations, reported confidence intervals, and worked with stability-assessed LST_cci climate data records. That study did not construct a trend by using only two early Terra years followed by an Aqua series. It also emphasized the desirability of records longer than 30 years for trend estimation. Similarly, regional MODIS trend work has used an Aqua-only, stability-assessed MYDCCI climate data record with a propagated uncertainty budget (https://doi.org/10.1080/01431161.2023.2240522). The methodological precedent therefore supports MODIS trend analysis when platform consistency and uncertainty are handled explicitly.
2. Physical distinction between LST, air temperature, and heat stress
Throughout the manuscript, radiometric skin temperature (LST), 2-m air temperature (T2m), and human heat stress are frequently treated as equivalent.
In Section 3.1, the authors assert that there is a constant offset of ~2.75 °C between daytime LST and maximum 2-m air temperature. The skin-to-air temperature difference is governed by surface energy balance partitioning and varies strongly across land covers: from near 0 °C in dense mountain forests to well over 15–20 °C over dry bare soils in the Bărăgan plain and impervious urban surfaces in Bucharest.
Moreover, fixed daytime LST thresholds (35, 40, 45, 50, 55 °C) are surface radiometric skin values and should not be labelled as physiological thresholds or air-temperature heatwaves. The authors also dismiss clear-sky sampling bias as having "negligible structural impact" (L. 182–185), yet acknowledge a mean cloud cover of 28.5% on warning days.
3. Structure and weighting of the Heat Vulnerability Index (HVI)
The index formulation (Section 2.3) presents several conceptual and mathematical issues:
4. Circular validation
In L. 701–706, the authors claim indirect validation of the HVI because high-HVI counties correlate with higher LST, more warning days, and larger elderly populations. This is a circular argument: showing that an index correlates with its own input variables, not independent empirical validation. To genuinely validate the index, it must be compared against external public health data (such as heat-related excess mortality or emergency medical calls), or else clearly framed as an unvalidated spatial screening index.
For all reported slopes, the manuscript must state the estimator, temporal unit, sample size, confidence interval, p-value, and handling of spatial multiple testing. Fig. 5 reports spatial percentages warming but does not identify statistically significant pixels. With only 26 summers, endpoint sensitivity and interannual variability are substantial. A simple pixelwise ordinary least-squares slope is insufficient.
Use a defensible approach such as Sen's slope with a modified Mann-Kendall test for autocorrelated series, or an equivalently justified model. Report field significance or false-discovery-rate control, trend confidence intervals, and sensitivity to influential years and the platform transition. The results should also be compared quantitatively with the longer air-temperature record rather than described as confirmation.
6. Missing regional literature
The manuscript claims to present the first national-scale heat vulnerability analysis for Romania, but overlooks foundational Romanian urban climate and satellite studies. In particular, the following works are highly relevant and must be integrated into the discussion:
- Cheval et al. (2022), MODIS-based climatology of the Surface Urban Heat Island at country scale (Romania), Urban Climate, 41, 101056, https://doi.org/10.1016/j.uclim.2021.101056. This earlier country-scale Romanian MODIS study addresses LST-air-temperature relationships, clear-sky limitations, and spatial controls.
- Cheval et al. (2023), A scale assessment of the heat hazard-risk in urban areas, Building and Environment, 229, 109892, https://doi.org/10.1016/j.buildenv.2022.109892. This is especially close: it combines MODIS LST, population density, and urban fabric across 77 Romanian cities.
- Mocanu et al. (2021), Human Health Vulnerability to Summer Heat Extremes in Romanian-Bulgarian Cross-Border Area, Natural Hazards Review, 22, https://doi.org/10.1061/(ASCE)NH.1527-6996.0000439. This develops a regional composite heat-health vulnerability index using exposure, sensitivity, and adaptive-capacity information.
- Cheval and Dumitrescu (2015), The summer surface urban heat island of Bucharest (Romania) retrieved from MODIS images, Theoretical and Applied Climatology, 121, 631-640, https://doi.org/10.1007/s00704-014-1250-8.
- Herbel et al. (2018), The impact of heat waves on surface urban heat island and local economy in Cluj-Napoca city, Romania, Theoretical and Applied Climatology, 133, 681-695, https://doi.org/10.1007/s00704-017-2196-4.
- Grigoras and Uritescu (2019), Land Use/Land Cover changes dynamics and their effects on Surface Urban Heat Island in Bucharest, Romania, International Journal of Applied Earth Observation and Geoinformation, 80, 115-126, https://doi.org/10.1016/j.jag.2019.03.009.
- Scripca, A.-S., Acquaotta, F., Croitoru, A.-E., and Fratianni, S. (2022), The impact of extreme temperatures on human mortality in the most populated cities of Romania, International Journal of Biometeorology, 66(1), 189-199, https://doi.org/10.1007/s00484-021-02206-w.
- Chitu, Z., Bojariu, R., Velea, L., and Van Schaeybroeck, B. (2023), Large sex differences in vulnerability to circulatory-system disease under current and future climate in Bucharest and its rural surroundings, Environmental Research, 234, 116531, https://doi.org/10.1016/j.envres.2023.116531.
- Zoran et al. (2026), Remote Sensing Monitoring of Summer Heat Waves-Urban Vegetation Interaction in Bucharest Metropolis, Atmosphere, 17, 109, https://doi.org/10.3390/atmos17010109. This also examines a long MODIS-era record through 2024.
Specific and Minor Comments
- L. 153–155: The statement that LST measures the "atmospheric greenhouse effect" is physically incorrect. Please rephrase in terms of radiometric surface temperature and thermal infrared emissions.
- Figure 1: This overview figure is presented before the Data and Methods section without explaining the underlying data source, spatial aggregation, or statistical regression.
- Figure 4: Please clearly specify the sample unit for n=2343 (are these station-days, county-days, or pixel-days?).
- Figures 3, 5, 6, 10, 11: The text and label sizes across multi-panel figures are much too small to read comfortably. Please enlarge all axis titles, tick labels, and legends.
- L. 376–377: The definition of heatwave days in 2024 (56 days) should be reconciled with the monthly sum of warning days (15 + 20 + 22 = 57 days).
- Wording / Style: Please replace informal or promotional phrases such as "empirical blueprint", "actionable diagnostic tool", and "transcends theoretical mapping" with sober scientific descriptions.
- Data and Code Availability: In line with Copernicus data policies, please deposit the processing scripts (GEE/Python/R) and processed figure-generating datasets in a permanent FAIR repository (e.g., Zenodo).
- The Conclusion is repetitive and substantially overstates operational and public-health applicability. Shorten it and include limitations and uncertainty alongside the results.
Minimum requirements for a credible new submission
1. Quantify the potential Terra-Aqua/observation-time discontinuity and incomplete 2002 season using Terra-only, Aqua-only, common-period, exclusion, and breakpoint tests; rebuild or restrict the series if the results show material bias.
2. Fully document QA, cloud/missing-data handling, aggregation, spatial processing, and trend inference.
3. Separate LST exceedances, air-temperature heatwaves, and operational warnings.
4. Restrict warning-trend conclusions to a homogeneous period or provide validated harmonization.
5. Redesign and rename the HVI/risk index; justify indicators, weights, normalization, and classes; quantify sensitivity and uncertainty.
6. Validate against independent health or impact outcomes, or explicitly present the output as an unvalidated screening index.
7. Archive code and figure-reproduction data; provide reviewer access to restricted inputs.
8. Reconstruct the literature review and verify every bibliographic record.