the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 04 Jan 2026
Forest fire risk in Croatia under a changing climate
Abstract. Forest fires in Croatia inflict substantial economic and ecological damage and frequently pose a threat to infrastructure and human lives. The southern part of the Croatian Adriatic, belonging to the Mediterranean basin, is the most severely affected region. To evaluate fire risk, the Canadian Fire Weather system was applied, and indices based on Fire Weather Index (FWI) – Seasonal Severity Rating (SSR), the number of days with FWI > 30 (FWI30), the 90-th percentile of FWI (FWIp90), and Length of Fire Season (LOFS) were derived. This study investigates the extent to which climate change has influenced the variability of latter indices across Croatia during June-September season. The analysis covers the period 1961–2020, revealing upward trends and predominantly positive anomalies in the evaluated indices. The most favourable fire weather conditions occur in the southern part of the Croatian Adriatic, which also exhibits the strongest increasing trends in SSR and FWI30. Although the continental parts of Croatia have historically been less susceptible to wildfires, the observed trends in the analysed indices suggest that conditions conducive to ignition and spread of wildfires are gradually emerging in these areas as well.
- Preprint
(2710 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 08 May 2026)
-
RC1: 'Comment on egusphere-2025-6439', Anonymous Referee #1, 15 Mar 2026
reply
-
AC1: 'Reply on RC1', Mislav Anić, 03 Apr 2026
reply
We would like to thank Referee #1 for her/his constructive suggestions and comments which have helped to improve the manuscript. We highly appreciate the time and effort invested in the review.
General comments
Anić et al. developed a 1-km dataset for fire-weather and related meteorological variables and indices based on station observations and digital elevation model (DEM) data. The results could be of interest to the local scientific community. My main concern is regarding the validity of the developed 1-km dataset and its added value.
Specifically, this dataset is essentially using 35 meteorological stations (24 Croatian + 11 Slovenian), 18 climatological stations, 105 rain gauges, interpolated to a 1 km grid by considering the elevation, distance from the sea, latitude, and longitude. Does the created dataset actually contain 1-km information? The network density is definitely much coarser than 1 km. On the other hand, even if information is added by elevation, this would mainly be relevant for temperature, while the relationship between elevation and precipitation/wind/cloud cover/sunshine duration/relative humidity is more complicated. It would be more convincing if the authors could add other sources of high-resolution data for climate and meteorological variables, such as satellite data. This is a common practice used for generating high-resolution observation datasets, such as for E-OBS and MSWEP.
The question of validity set aside, what added value has the developed high-resolution dataset provided? What part of the conclusion is not what we’ve already known based on other observation, satellite, reanalysis, or model datasets?
We agree with Referee #1 that the station density is relatively coarse. If we had analysed the period 1981–2020, data from more stations would have been available, however during the period 1981–2020, a warmer climate regime is already strongly expressed. Since the aim of this study was to analyse trends, we opted for the longer period 1961–2020, as a longer time series provides a more reliable signal of climate change and reduces the influence of natural variability. We are aware that the success of regression kriging depends, among other factors, on the density and spatial distribution of stations used for interpolation. To address the aim of the study, two approaches were at our disposal: 1) to calculate the FWI using spatially interpolated daily meteorological values, or 2) to calculate FWI from the station data, then aggregate them seasonally, and subsequently calculate climate normal, and afterwards perform spatial interpolation of seasonal values and climate normal. Considering that daily values exhibit substantially higher variability than seasonal values, we decided to use approach 2. The fact that the interpolation was based on indices rather than the original meteorological data actually contributed to the robustness of the estimates. We agree that this methodology has certain limitations, nevertheless we consider that it still provides valuable insight into changes in fire danger in Croatia during the period 1961–2020. Our objective was not to develop and validate high-resolution meteorological dataset, but rather to assess the extent to which climate change has influenced variations in FWI-based indices across the territory of Croatia. We also agree that the use of satellite data could improve the spatial representation and accuracy of the estimates. However, for the period considered in our study, particularly the earlier decades (1961–1990) satellite-based data are either unavailable or limited in spatial resolution. Nevertheless, we acknowledge the potential of satellite data and plan to include them in future work, especially for higher-resolution analyses and more recent periods. The results obtained in this study are consistent with previous studies conducted at the Mediterranean or European scale, although those studies were performed at rather coarser spatial resolutions than 1 km. Croatia is characterized with complex topography, and in some areas elevation changes rapidly, meaning that even a 1 km resolution may still be relatively coarse in certain locations. Croatia is characterized by highly complex terrain, with pronounced elevation gradients over short distances, as well as the coexistence of Mediterranean, mountainous, and continental climate regimes within a relatively small area. In such cases, coarse-resolution datasets tend to smooth spatial gradients and may therefore underestimate local extremes and variability. With the presented approach, we aimed to assess the evolution of fire danger at a finer spatial scale, explicitly accounting for the complex interplay between climate and topography in the study area. Although the results are based on spatially interpolated indices, the increased spatial resolution allows for the identification of localized patterns and changes in fire danger that are not resolved in coarser datasets. This is particularly important for areas where fire danger has historically been low but is now increased and existing fire protection system may need to be reconsidered. The results therefore provide additional, spatially detailed information for Croatia, while also contributing to a broader Mediterranean context, particularly in transitional zones between climatic regimes and in topographically complex areas.
In in the revised manuscript within the discussion, we will add a paragraph in which we critically reflect on the limitations of the method used for spatial interpolation.
Specific comments
Lines 39-40: What are the causes of wildfire ignition in JJAS (the second peak)?
The majority of wildfires in Croatia are associated with human activity (Mataković et al., 2024). During the JJAS season, there is a substantial increase in the number of people staying in Croatia, especially in its coastal regions. Previous research showed that the most common causes of wildfire ignition are arson and negligence by both domestic and foreign visitors and tourists, while natural causes such as lightning account for a smaller proportion of cases.
In a revised version of the manuscript, we will add this explanation in chapter 1 Introduction.
Mataković, H., Beljan, K., and Posavec, S.: Perception of the causes and consequences of forest fires in the Republic of Croatia (in Croatian), Šumar. List, 148(7-8), 327-340. https://doi.org/10.31298/sl.148.7-8.1, 2024
Lines 126-129: Please provide references for the definition of daily severity rating. Why is it defined this way? How are the values for the coefficient and exponent chosen? What is the physical meaning of this quantity? How is its physical meaning different from FWI?
To assess the fire risk, the Canadian Fire Weather Index (FWI) system was applied in this study. The calculations followed the methodology described by Van Wagner (1987). According to Van Wagner, the Daily Severity Rating (DSR) was introduced as an optional component of the FWI system and is computed directly from the FWI index using the equation DSR = 0.0272 · FWI¹·⁷⁷. The FWI represents the intensity of a spreading fire as the energy output rate per unit length of the fire front and is based on the rate of fire spread and the amount of available fuel. The DSR weights the FWI as its values increase, in a way that is considered to reflect the difficulty of fire control in more direct proportion to the amount of work required to suppress a fire. Furthermore, according to Van Wagner, the FWI itself is not considered suitable for averaging and should be used only as a daily value. Any averaging, whether spatially across multiple stations on a given day or temporally at a single station over a certain period, is better performed using the DSR. For this reason, seasonal averaging of the FWI and spatial interpolation of the FWI index were not performed. Instead, the DSR was seasonally averaged, which allowed the calculation of seasonal severities (SSR), and spatial interpolation was subsequently performed for the SSR. Furthermore, the severity rating was originally introduced by Williams (1959) to provide a measure of control difficulty in terms of the former fire danger index. The exponent 1.77 was obtained through empirical fitting during the development of the FWI system, and the power transformation of the FWI index was selected because it reproduces well the nonlinear increase in fire severity. The parameter 0.0272 is a calibration factor with the purpose to rescale the DSR into a practical numerical range.
We will improve section 3.2 FWI calculation with more details about daily severity rating in the revised manuscript.
Williams, D. E.: Fire season severity rating, Forestry Branch, Department of Northern Affairs and National Resources, Forest Research Division, Technical Note 73, 13 pp., Ottawa, Canada, 1959.
Van Wagner, C. E.: Development and structure of a Canadian forest fire weather index system, Forestry Technical Report 35, Canadian Forestry Service, Ottawa, 1987.
Line 145: Is the leave-one-out cross-validation for leaving out one station or one or more years? Please provide more details on what (which stations/years) are used for model training, and what is left for validation.
Leave-one-out cross-validation (LOOCV) was performed by omitting one station at a time, after which a kriging prediction was made for the omitted station based on the defined model. The residual was then calculated as the difference between the observed and predicted values. This procedure was repeated for each station used in the spatial interpolation. LOOCV was not carried out by omitting individual years or seasons; instead, a regression kriging prediction was performed for each season, and LOOCV was applied to each of these predictions in the manner described above. Thus, in the case of predicting seasonal values, we obtained 60 values of the coefficient of determination (R²), mean error (ME), and normalized root mean square error (RMSEr). From these values, the mean, minimum, and maximum were calculated and are presented in Table 1 in the row seasonal values. In the case of climate normals, a single LOOCV was performed, and the corresponding R², ME, and RMSEr are shown in the same table in the row normals.
We will improve section 3.3. Spatial interpolation with more details about LOOCV in the revised manuscript.
Lines 296-298: Why does a larger number of stations lead to higher error? Wouldn’t more information be added? Could this larger error be simply due to the larger variability at the seasonal timescale?
A larger number of stations (157) used for spatial interpolation were employed exclusively in the case of the number of dry days (Ndry) and dry spell durations (Ddry), as additional rain gauge stations measuring only precipitation were available. Because precipitation exhibits higher spatial variability than other variables, the correlation with predictors is weaker compared to, for example, maximum air temperature. Therefore, estimation error is larger compared to the estimation of some other observed variables, especially maximum air temperature. Nevertheless, the increased station density provided additional information, improving the interpolation of the number of dry days and dry spell duration. In the case of interpolating seasonal values of Ndry and Ddry, the error was higher than for their climatological normals, due to greater variability at the seasonal time scale.
We will add this information in section 4.3 Accuracy of spatial prediction in the revised manuscript.
Lines 332-334: The negative correlation between temperature and precipitation is merely a statistical observation. I don’t think it is sufficient to be the driver of the positive land-atmospheric feedback. Please provide physical explanations for the positive feedback mechanism.
Thank you for this comment. We agree that the negative correlation between temperature and precipitation alone does not provide a physical explanation for the positive land–atmosphere feedback. Specifically, soil moisture deficits limit evapotranspiration and latent heat flux, while increasing sensible heat flux and near-surface air temperature. At the same time, reduced evapotranspiration may suppress precipitation, further reinforcing dry conditions. Conversely, enhanced atmospheric demand during heatwaves can accelerate soil moisture depletion, resulting in self-amplifying feedback between temperature extremes and drought (Hao et al., 2022). We plan therefore clarify the underlying physical mechanisms within Discussion in the revised manuscript.
Hao, Z., Hao, F., Xia, Y., et al.: Compound droughts and hot extremes: Characteristics, drivers, changes, and impacts, Earth-Sci. Rev., 235, 104241, https://doi.org/10.1016/j.earscirev.2022.104241, 2022.
Technical corrections
Line 278: “the FWI30 and FWIp90 trend results” → “the FWI30 trend and FWIp90 anomalies”
Corrected.
Line 360: “…, could potentially…” → “…, which could potentially…”
Corrected.
Citation: https://doi.org/10.5194/egusphere-2025-6439-AC1 -
RC2: 'Reply on AC1', Anonymous Referee #1, 10 Apr 2026
reply
The approach 2 (calculate FWI at station level, aggregate temporally (or is DSR the variable that is aggregated?), and then interpolate) is fine. And I have no doubt that increased spatial resolution would benefit the fire danger evaluation and management in Croatia. The issue I raised is whether the interpolated data obtained by the authors actually contains the amount of information they claimed. I see in the manuscript that the authors still aim to “assess … at 1 km spatial resolution.”
More specifically:
Does providing a 1-km DEM make the interpolated FWI actually contain 1-km information?
Does FWI variation follow the terrain as much as temperature does? I don’t think DEM alone is sufficient to explain the spatial variability of temperature, not to mention FWI, which also depends on relative humidity, wind speed, and precipitation.
Does the terrain-adjusted FWI at a grid agree with the actual FWI that would have been measured there? How did you verify this?Citation: https://doi.org/10.5194/egusphere-2025-6439-RC2 -
AC2: 'Reply on RC2', Mislav Anić, 15 Apr 2026
reply
The approach 2 (calculate FWI at station level, aggregate temporally (or is DSR the variable that is aggregated?), and then interpolate) is fine. And I have no doubt that increased spatial resolution would benefit the fire danger evaluation and management in Croatia. The issue I raised is whether the interpolated data obtained by the authors actually contains the amount of information they claimed. I see in the manuscript that the authors still aim to “assess … at 1 km spatial resolution.”
We are grateful to the Referee #1 for pointing at this issue. The Referee #1 is correct that the interpolated values do not carry the amount of information as a reader might assume based on the wording in our aims from the original submission.
We will amend the wording of the Aim (L62-63) with the following sentence to assure clarity:
“The aim of this study is to assess the trends of wildfire risk by analysing spatial and temporal variability in climate and FWI based indices in Croatia during the period 1961–2020, with indices derived from station observations and interpolated to a 1 km grid for spatial representation.”
We first calculated the FWI index following the methodology described in Van Wagner (1987), after which we derived several FWI-based indices, including DSR. FWI-based indices were then temporally aggregated (summed in the case of the number of days with FWI > 30) and mapped onto an interpolation grid with a horizontal resolution of 1 km using the regression kriging method.
Van Wagner, C. E.: Development and structure of a Canadian forest fire weather index system, Forestry Technical Report 35, Canadian Forestry Service, Ottawa, 1987.
More specifically:
Does providing a 1-km DEM make the interpolated FWI actually contain 1-km information?
We propose the amendment to the wording of the aims, which clarifies this issue (see above). In addition, we agree that the available number of stations that was used for spatial interpolation with regression kriging in this study is relatively small, which represents a limitation of this approach. However, since our objective was to assess trends in FWI-based indices over Croatia, we opted to analyse a longer period (1961–2020), for which data from a smaller number of stations were available. Nevertheless, the evaluation of model performance indicated that the predictions are reasonably reliable.
In the revised manuscript, we will add the following text to the discussion, providing a critical assessment of the prediction:
„In this study, regression kriging was applied to interpolate FWI-based indices calculated from meteorological station data and express them on a 1 km grid, providing a spatially continuous representation of regional patterns to assess changes in fire danger over Croatia during the period 1961–2020. The main limitation of this approach is the relatively small number of stations available for interpolation in this study. Despite this, the model performance evaluation based on LOOCV indicated that the predictions are reliable and provide valuable insight into changes in fire danger across Croatia over the study period. The prediction was more successful for the interpolation of climatological normals than for seasonal values, which is expected given that seasonal values exhibit higher variability compared to normals and therefore show weaker correlations with the predictors. The use of a larger number of stations would undoubtedly provide more information and improve the accuracy of the model prediction.“
Does FWI variation follow the terrain as much as temperature does? I don’t think DEM alone is sufficient to explain the spatial variability of temperature, not to mention FWI, which also depends on relative humidity, wind speed, and precipitation.
We agree with the Reviewer that interpolation of FWI is challenging, precisely because FWI depends on multiple meteorological variables and characteristic of available biomass (i.e. fuel), which may not be directly related to elevation or other terrain characteristics. Currently there are no reliable wind and precipitation data at 1 km resolution for the period of our investigation. Creating such high-resolution maps at 1 km would be highly challenging and it was outside the scope of this manuscript. However, we are aware of the limitations, and will add in the revised manuscript a text which reflects on the limitations of the FWI based indices interpolation based on the method we used.
Considering these limitations, we decided to apply the conservative approach and use standard predictors commonly applied in climatology for spatial interpolation, namely elevation, latitude, longitude, and distance from the sea (e.g., Hiebl et al., 2009; Perčec Tadić et al., 2010, 2022; Crespi et al., 2019; Omazić et al., 2024). Leave-one-out cross-validation showed that the prediction was most successful for maximum air temperature (Tmax), where (in case of climate normal) the model explained 88% of the spatial variability (R² = 0.88). For the FWI-based indices, the prediction performance was lower: the model explained 64% of variability for FWIp90, 67% for FWI30, 79% for SSR, and 82% for LOFS. As expected, prediction performance was better for climatological normals, since seasonal values exhibit higher variability than normals, resulting in weaker correlations between seasonal values and the predictors. Consequently, performance metrics (R², ME, RMSEr) were lower for the interpolation of seasonal values.
Perčec Tadić, M.: Gridded Croatian climatology for 1961–1990, Theor. Appl. Climatol., 102, 87–103, https://doi.org/10.1007/s00704-009-0237-3, 2010.
Perčec Tadić, M., Pasarić, Z., and Guijarro, J. A.: Croatian high-resolution monthly gridded dataset of homogenised surface air temperature, Theor. Appl. Climatol., 151, 227–251, https://doi.org/10.1007/s00704-022-04241-y, 2023.
Hiebl, J., Auer, I., Böhm, R., Schöner, W., Maugeri, M., Lentini, G., Spinoni, J., Brunetti, M., Nanni, T., Perčec Tadić, M., Bihari, Z., Dolinar, M., and Müller-Westermeier, G.: A high-resolution 1961–1990 monthly temperature climatology for the greater Alpine region, Meteorol. Z., 18, 507–530, https://doi.org/10.1127/0941-2948/2009/0403, 2009.
Crespi, A., Lussana, C., Brunetti, M., Dobler, A., Maugeri, M. And Tveito, O. E: High-resolution monthly precipitation climatologies over Norway (1981–2010): Joining numerical model data sets and in situ observations, Int. J. Climatol., 39, 2057–2070, DOI: 10.1002/joc.5933, 2019.
Omazić, B., Anić, M., Telišman Prtenjak, M., Kvakić, M., and Blašković, L.: Analysis of different existing measurement-based methods and a new approach for frost probability detection, Agric. For. Meteorol., 347, 109898, https://doi.org/10.1016/j.agrformet.2024.109898 , 2024.
Does the terrain-adjusted FWI at a grid agree with the actual FWI that would have been measured there? How did you verify this?
To assess model performance, we performed leave-one-out cross-validation (LOOCV). By omitting one station at a time and predicting its value based on the remaining stations, LOOCV provided an estimate of model performance at locations without observations. In this way, we effectively conducted a validation at locations without measurements.
As already mentioned in our previous replies, we will expand the discussion on the limitations of the applied interpolation method.
Citation: https://doi.org/10.5194/egusphere-2025-6439-AC2 -
RC3: 'Reply on AC2', Anonymous Referee #1, 16 Apr 2026
reply
The authors have sufficiently addressed my comments. I have no further comments.
Citation: https://doi.org/10.5194/egusphere-2025-6439-RC3
-
RC3: 'Reply on AC2', Anonymous Referee #1, 16 Apr 2026
reply
-
AC2: 'Reply on RC2', Mislav Anić, 15 Apr 2026
reply
-
RC2: 'Reply on AC1', Anonymous Referee #1, 10 Apr 2026
reply
-
AC1: 'Reply on RC1', Mislav Anić, 03 Apr 2026
reply
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 303 | 179 | 23 | 505 | 14 | 17 |
- HTML: 303
- PDF: 179
- XML: 23
- Total: 505
- BibTeX: 14
- EndNote: 17
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
General comments
Anić et al. developed a 1-km dataset for fire-weather and related meteorological variables and indices based on station observations and digital elevation model (DEM) data. The results could be of interest to the local scientific community. My main concern is regarding the validity of the developed 1-km dataset and its added value.
Specifically, this dataset is essentially using 35 meteorological stations (24 Croatian + 11 Slovenian), 18 climatological stations, 105 rain gauges, interpolated to a 1 km grid by considering the elevation, distance from the sea, latitude, and longitude. Does the created dataset actually contain 1-km information? The network density is definitely much coarser than 1 km. On the other hand, even if information is added by elevation, this would mainly be relevant for temperature, while the relationship between elevation and precipitation/wind/cloud cover/sunshine duration/relative humidity is more complicated. It would be more convincing if the authors could add other sources of high-resolution data for climate and meteorological variables, such as satellite data. This is a common practice used for generating high-resolution observation datasets, such as for E-OBS and MSWEP.
The question of validity set aside, what added value has the developed high-resolution dataset provided? What part of the conclusion is not what we’ve already known based on other observation, satellite, reanalysis, or model datasets?
Specific comments
Lines 39-40: What are the causes of wildfire ignition in JJAS (the second peak)?
Lines 126-129: Please provide references for the definition of daily severity rating. Why is it defined this way? How are the values for the coefficient and exponent chosen? What is the physical meaning of this quantity? How is its physical meaning different from FWI?
Line 145: Is the leave-one-out cross-validation for leaving out one station or one or more years? Please provide more details on what (which stations/years) are used for model training, and what is left for validation.
Lines 296-298: Why does a larger number of stations lead to higher error? Wouldn’t more information be added? Could this larger error be simply due to the larger variability at the seasonal timescale?
Lines 332-334: The negative correlation between temperature and precipitation is merely a statistical observation. I don’t think it is sufficient to be the driver of the positive land-atmospheric feedback. Please provide physical explanations for the positive feedback mechanism.
Technical corrections
Line 278: “the FWI30 and FWIp90 trend results” → “the FWI30 trend and FWIp90 anomalies”
Line 360: “…, could potentially…” → “…, which could potentially…”