| 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.
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.