Machine learning assessment of climate-driven variability in European forest fire burnt areas

RESEARCH PAPER https://doi.org/10.1071/xxxx Machine learning assessment of climate-driven variability in European forest fire burnt areas Emil CiencialaA,B,* , Milan FischerB,C , Lucie KudlᡠckováB,C , Markéta Podˇ ebradskáB,C , Jan BalekB,C , Radka MaškováA, Petr Štˇ epánekB, Jana BeranováA,and Miroslav TrnkaB,C ABSTRACT Background. Despite intensifying fire-conducive weather, European burnt forest areas show contradictory trends, suggesting climate impacts may be moderated by fire management systems. Aims. We tested whether climate variables remain the primary drivers of burnt forest area variability across Europe, and whether machine learning models can robustly predict fire outcomes despite three decades of management improvements. Methods. We analyzed forest fire data from 12 European countries (1990–2022) from UNFCCC reports, developing ensemble machine learning models using fire weather indices, dead-fuel moisture, soil moisture, and country as a categorical variable to capture management differences. Key results. The ensemble model explained 84% and 80% of variance in calibration and verification datasets. Despite increased fire-conducive weather in six countries, 10 of 12 showed stable or declining burnt areas, with Poland and Spain declining significantly. Czech Republic projections under SSP2-4.5 scenarios indicate burnt area may increase 1.5–2.0-fold by 2085. Conclusions. Climate-based models remain robust predictors of fire outcomes. European fire management has successfully decoupled outcomes from climate pressure, but future climate change may challenge this effectiveness. Implications. Adaptive management strategies and proactive knowledge transfer from Mediterranean to temperate regions are crucial in the face of progressing climate change. Keywords: climate change, fire risk, fire weather, forest fires, machine learning. Introduction The changing climate and increasing frequency and severity of1 forest fires in recent years (Jones et al.,2024;Zhang et al.,2024)2 have highlighted the importance of fire prevention and fire mit-3 igation strategies. For these strategies, information on the likely4 development of fire risk and the extent of forest fires associated5 with a changing climate is essential. However, critical ques-6 tions remain: How predictable are future burnt areas based on7 climate drivers? And can European fire management systems8 sustain their effectiveness under continued climate change?9 Addressing these questions requires understanding both histor-10 A. IFER - Institute of Forest Ecosystem Research, ˇ Cs. armády 655, Jílové u Prahy, 254 01, Czech Republic, https://orcid.org/0000-00021254-4254; B. Global Change Research Institute CAS, Bˇ elidla 986/4a, Brno, 603 00, Czech Republic, https://orcid.org/0000-0002-7841-9317 https://orcid.org/0000-0002-3921-1601 https://orcid.org/0000-0002-31214904 https://orcid.org/0000-0002-2753-3777 https://orcid.org/00000001-8956-5590 https://orcid.org/0000-0003-4727-8379 ; C. Department of Agrosystems and Bioclimatology, Faculty of Agronomy, Mendel University in Brno, Zemˇ edˇ elská 1, Brno, 613 00, Czech Republic, https://orcid.org/0000-0002-7841-9317 https://orcid.org/0000-00023921-1601 https://orcid.org/0000-0002-3121-4904 https://orcid.org/00000002-2753-3777 https://orcid.org/0000-0003-4727-8379 *Corresponding author. Email: [email protected] ical fire management success and future fire risk under climate 11 scenarios. 12 Improved fire monitoring and early interventions in the 13 event of fire events help offset the effects of a changing cli14 mate and generally mitigate fire risk. In the face of climate 15 change, fire activity is intensifying in many regions, particularly 16 where hot and dry conditions are becoming more frequent (for 17 example, Abatzoglou et al. (2025)). Understanding how climate 18 change translates into actual burnt area is crucial for evaluating 19 the true extent of wildfire risk. 20 However, the relationship between climate drivers and burnt 21 area is not universally linear or straightforward (Burton et al.,22 2024), and a critical paradox has emerged in European fire 23 regimes. While numerous studies have reported increasing 24 trends in fire-conducive weather - including higher temper25 atures, reduced precipitation, and increased vapor pressure 26 deficit - actual burnt areas do not always follow the same tra27 jectory. Remarkably, several Mediterranean countries, among 28 Europe’s most fire-prone regions, have demonstrated stable or 29 even decreasing trends in burnt forest area despite increasingly 30 arid climatic conditions (Turco et al.,2016;Jones et al.,2022). 31 This apparent paradox reflects the growing effectiveness of 32 national and regional fire management systems, which include 33 early warning networks, suppression infrastructure, prescribed 34 burning, and land use planning. As a result, while burnt areas 35 Cienciala et al. International Journal of Wildland Fire typically increase under warming conditions, they may decline36 in regions that implement coordinated and sustained fire pre-37 vention strategies (Schulte et al.,2017). This disconnect between38 fire-prone weather and observed fire outcomes highlights the39 need for a modeling framework that integrates both climatic40 drivers and the effectiveness of national fire management sys-41 tems when assessing and projecting wildfire risk across Europe.42 Data-driven methods, particularly machine learning ensem-43 ble approaches, offer considerable promise for this task: they44 can quantify the relative importance of climate variables while45 capturing complex, nonlinear interactions between fire drivers46 and outcomes, and they provide operational tools for fire man-47 agement agencies to project future fire risk under climate sce-48 narios.49 Moreover, the role of climate as a fundamental driver of fire50 dynamics remains evident, particularly through the fire weather51 index (FWI, Wagner (1987)), vapor pressure deficit (VPD), and52 dead fuel moisture content (DFMC, de Dios et al. (2015)), etc.53 These variables are widely used in fire risk assessment and have54 shown high predictive power for fire occurrence and spread55 (Venäläinen et al.,2014;Resco de Dios et al.,2021;Kang et al.,56 2023;Wang et al.,2023;Nagavciuc et al.,2024). These variables57 and others are further complemented by remote-sensing obser-58 vations, which provide valuable spatial and temporal insights59 into actual fire activity across forested landscapes.60 There is a range of remote-sensing approaches to estimate61 burnt forest areas from properties applicable on a fine time62 scale to assess the spread and development of forest fires locally63 (Siljander,2009;Pessôa et al.,2020). While remote sensing tech-64 niques provide critical data for mapping burnt areas and track-65 ing fire events at high spatial and temporal resolutions, there66 is increasing demand for scalable statistical models capable of67 capturing long-term patterns and enabling comparative fire risk68 assessments at national to continental levels. Such approaches69 are particularly important for policy-relevant fire risk assess-70 ments in Europe, where strong gradients in climate, forest71 composition, and land-use governance require harmonized yet72 flexible modeling strategies.73 In response to this need, we propose a climate-driven ensem-74 ble modeling framework that incorporates fire-relevant bio-75 physical predictors and machine learning techniques to esti-76 mate burnt forest area across 12 European countries, analyzing77 fire patterns at a country to regional scales across the European78 gradient.79 The central hypothesis guiding this study is that climate-80 based driving forces have remained decisive for the develop-81 ment of burnt forest areas in Europe since 1990. We ask whether82 climate-based statistical models remain sufficiently robust in a83 changing technological and climatic landscape to inform prac-84 tical, forward-looking wildfire risk assessments. To explore this85 question, we analyze the observed forest fires in 12 European86 countries, using the burnt forest area information reported in87 the national GHG emission inventories to the United Nations88 Framework Convention on Climate Change (UNFCCC), com-89 bined with climate-based drivers and forest resource informa90 tion. 91 These specific objectives were defined as follows: i. To eval92 uate the robustness of climate-based statistical models for pre93 dicting burnt forest area across 12 European countries for the 94 period of 1990–2022; ii. To test the applicability of such mod95 els for projecting future fire risk expressed as the Burnt Forest 96 Area Index (BFAI) under evolving climate scenarios for the 97 remainder of the 21st century. 98 Methods Project domain and temporal scale 99 The analyses carried out for this study included 12 European 100 countries, including Austria (AUT), Bulgaria (BGR), the Czech 101 Republic (CZE), Spain (ESP), France (FRA), Greece (GRC), 102 Croatia (HRV), Italy (ITA), Poland (POL), Portugal (PRT), 103 Slovenia (SVN) and Sweden (SWE). This selection of countries 104 represents a South-North European gradient (Fig. 1). It omits 105 coastal regions and countries with limited spatial extent, which 106 interferes with the sample size. The assembled data on forest 107 fire areas included country-level data on burnt forest areas from 108 1990 to 2022 with annual resolution. Fig. 1. Analyzed domain of the study, including 12 European countries. 109 Source data on forest fires and forest 110 resources 111 Data on forest area (ha/year), forest fire area (ha/year) and 112 estimated emissions (CO2eq./year) attributed to forest fires 113 were taken from the Common Reporting Format (CRF) 114 2 www.publish.csiro.au/xx International Journal of Wildland Fire Fig. 2. Box plot (with outliers shown by red points) of the annual burnt areas for the 1990–2022 period of the 12 European countries, distinguishing the calibration, verification, and entire datasets. Note: A log scale is used on the y-axis. tables, which are a mandatory part of the national inven-115 tory submissions to UNFCCC (UNFCCC,2024), together with116 the National Inventory Report texts (https://www.unfccc.int;117 accessed 30th June 2024). We also checked the fire area data118 from the European Forest Fire Information System (European119 Commission,2024), which do not fully correspond to the120 UNFCCC data. Therefore, UNFCCC data were used exclusively121 in this study, as reported fires are attributed to forest land, and122 all the NIR reports are regularly verified during the UNFCCC123 review process. In addition to the burnt forest area, we also col-124 lected the extent of forest area per country from the national125 CRF tables. The main dependent variable was the annual burnt126 forest area, expressed as a fraction related to 1000 ha (kha) of127 the total forest area in the country and referenced as the BFAI128 (ha/kha). This dependent variable is shown in Fig. 2, both for129 theentiredatasetandfor specific subsets, as describedinSection130 2.3 below. Finally, to analyze the possible effects of forest tree131 species composition, we used a tree species map for Europe132 (Brus et al.,2012). This map provides gridded (1 ×1 km) infor-133 mation on 20 tree species groups in Europe. It was overlaid on134 the 12 analyzed European countries, and the share of aggregated135 broadleaved, coniferous, and specifically pines species groups136 was retrieved. This information served as additional indepen-137 dent data in addition to the climatic and climate-related indices138 described below.139 Forest fire-related indices 140 The selected indices presumably related to forest fire risk 141 included I) the mean FWI in the form of a dimensionless 142 index of potential forest fire intensity, calculated from daily 143 temperature, relative humidity, wind speed, and precipitation 144 ((Wagner,1987), FWI_DI; -); II) FWI (fire danger classes accord145 ing to Kudláčková et al. (2024), where 3 is moderate, 6 is 146 extreme) as the number of days with an FWI greater than 3, 147 4, 5 and 6 (FWI_3+, FWI_4+, FWI_5+, FWI_6; days); III) 148 mean of the 10-hour DFMC (according to Fire Weather Indices 149 WIKI (nd); DFMC10H; %); IV) number of days with dead fuel 150 moisture less than 10%, 8%, 6%, 4% and 2% (DFMC10H_u10, 151 DFMC10H_u8, DFMC10H_u6, DFMC10H_u4, DFMC10H_u2; 152 days); V) the vapor pressure deficit (VPD, kPa), computed from 153 mean daily air temperature and relative humidity using the 154 standard Clausius–Clapeyron relation; VI) available relative 155 water content (AWR) as the number of days with soil mois156 ture under 30% for the 0–10 cm and 0–40 cm soil profiles 157 (AWR-10_u30, AWR-40_u30; days); VII) number of days with 158 available water percentiles (AWPs) in the 0–40 cm soil profile 159 above 2, 3 and 4 (AWP0_40_S2+, where AWR < 10th percentile; 160 AWP0_40_S3+, where AWR < 5th percentile; AWP0_40_S4+, 161 where AWR < 2nd percentile; days), and the AWR represents 162 the drought intensity for a given day relative to ±10 days in the 163 1981–2020 period; and VIII) available cumulative water deficit 164 (AWD) in the 0–40 cm soil profile (AWD0_40_sum; mm), where 165 AWD is the soil moisture deficit compared with normal for a 166 given day relative to ±10 days in the 1981–2020 period and 167 3 Cienciala et al. International Journal of Wildland Fire is calculated as the difference between the AWR and the 50th168 percentile. The AWP, AWD and AWR were calculated accord-169 ing to Hlavinka et al. (2011) and Trnka et al. (2020). For each of170 the 18 independent variables under I) to VIII), aggregates from171 daily values were created for the entire year, and specific sub-172 periods were determined by months, namely March to October173 (MAMJJASO), April to September (AMJJAS), October to March174 (ONDJFM), and November to February (NDJF). This resulted175 in 90 independent variables, considered in the statistical anal-176 ysis, in addition to the shares of broadleaves (or coniferous)177 and pines that were treated as constants at the country level178 for the entire period. Finally, to account for the different (abi-179 otic) factors linked to specific countries, we considered country180 as a categorical predictor. For the analysis, the input data were181 split into calibration and verification datasets, representing 2/3182 and 1/3 of the input data, respectively. The calibration dataset183 included the first and second years of each 3-year, nonoverlap-184 ping block within the analyzed period from 1990 to 2022, while185 theverificationdataset consisted of every thirdyearwithin those186 blocks. The resulting values of the dependent variable BFAI for187 the calibration, verification, and entire datasets are shown in188 Fig. 2. The Kruskal–Wallis test confirmed no significant differ-189 ences in BFAI distributions between calibration and verification190 datasets across all countries (p > 0.05).191 Machine learning approach for estimating the192 BFAI193 Several statistical approaches were used to test whether the194 BFAI is predictable on the basis of fire weather variables and195 basic parameters, such as country and forest type. All 12 coun-196 tries spanning a large gradient of environmental conditions197 across Europe were used for this analysis. Owing to the nonnor-198 mal distribution of the BFAI across all countries, the data were199 log-transformed before all analyses.200 Feature selection and dimensionality reduction201 Given the high dimensionality of the predictor variables (cfr.202 ’Forest fire-related indices’ subsection), a multistep feature203 selection procedure was implemented.204 1. Correlation filtering was applied to remove multicollinear205 predictors. Variables with pairwise correlation coefficients206 exceeding 0.9 were examined, and in each highly correlated207 pair, the variables with weaker correlations to the response208 variable were discarded. This step left 24 out of the original209 94 predictors.210 2. Boruta selection (Kursa and Rudnicki,2010), an algorithm211 that identifies all relevant features by comparing their impor-212 tance to that of randomized shadow variables, was applied213 using 10-fold cross-validation. Predictors selected in at least214 five out of ten folds were retained. This step ensured robust215 and objective identification of 15 consistently important fea-216 tures.217 3. Recursive feature elimination was conducted using random 218 forests to rank variable importance (Kuhn and Johnson,219 2013). As a part of final feature engineering, we decided to 220 exclude variables of the same type if they were moderately 221 correlatedand/or hadsimilar importance (e.g., VPD and FWI 222 for different aggregation periods). As a consequence, nine 223 predictors were kept for further analysis. 224 Model specification and preprocessing strategies 225 A total of 12 regression algorithms were configured using the 226 tidymodels framework (Kuhn and Wickham,2020;Kuhn and 227 Silge,2022). These included: 228 •Linear models: Generalized Linear Model (GLM) and Elastic 229 Net. 230 •Flexible regression: Multivariate Adaptive Regression Splines 231 (MARS). 232 •Kernel-based methods: Support vector machines with linear, 233 polynomial, and radial kernels, 234 •Tree-based models: Random Forest, Bagged Trees, 235 Classification and Regression Trees (CART), and Cubist. 236 •Boosting and additive models: XGBoost and Bayesian 237 Additive Regression Trees (BART). 238 •Artificial neural networks: Multilayer Perceptron (MLP). 239 Each algorithm was combined with one of three preprocess240 ing strategies: (i) no preprocessing, (ii) normalization of numer241 ical predictors, or (iii) polynomial expansion and interaction 242 terms among numeric features. 243 Country was included as a categorical predictor using like244 lihood encoding based on a logistic regression fitted to the 245 response variable. This approach preserves ordinal relation246 ships between the country and the outcome while avoiding the 247 increase in dimensionality of one-hot encoding. 248 Workflows combining preprocessing and model specifica249 tion were constructed and stored for tuning. 250 Model tuning and evaluation 251 Two hyperparameter tuning strategies were employed to ensure 252 both robustness and computational efficiency. First, a grid 253 search method with a fixed grid of 50 parameter combinations 254 per model was used. Second, a racing (ANOVA) method with an 255 early-stopping algorithm was used to evaluate up to 100 config256 urations per model. The inclusions of both methods allowed for 257 cross-verification of optimal model settings and ensured that the 258 selected configurations were not artifacts of any single tuning 259 strategy. Tuning was performed using 10-fold cross-validation 260 on the calibration dataset. Models were ranked using the root 261 mean square error (RMSE), and the best-performing config262 uration for each model was selected. The performance was 263 evaluated separately on the calibration and verification datasets. 264 Predictions were generated using the finalized models, and per265 formance metrics (RMSE and R2) were computed to assess 266 model accuracy and generalizability. 267 4 www.publish.csiro.au/xx International Journal of Wildland Fire Ensemble model construction268 To improve predictive robustness, an ensemble model was con-269 structed using stacked generalization. Predictions from the top-270 performing individual models were combined via penalized271 linear regression (stacking). The stacking coefficients were opti-272 mized using cross-validation, and the ensemble process was273 repeated20 timeswithvaryingseedstoensurestability.Thecon-274 figuration yielding the lowest calibration RMSE was selected as275 the final ensemble model.276 Because penalized regression does not force just a single con-277 figuration per algorithm, multiple tuned variants of the same278 model type might be retained in the final stack if they provide279 complementary predictive information.280 Projection of the BFAI281 For a detailed examination of the 21st-century climate under282 socioeconomic development conditions as represented by the283 SSP2-4.5 scenario, we employed the most recent CMIP6 global284 climate model (GCM) projections (Eyring et al.,2016;O’Neill285 et al.,2016). Daily data from 18 simulations across 17 GCMs286 were sourced from the Earth System Grid Federation (Cinquini287 et al.,2014), with most simulations at a nominal 100 km resolu-288 tion, one model at 50 km, and a subset at a coarser 250 km reso-289 lution. For the Czech Republic, CMIP6 GCMs were validated290 across various meteorological variables, including daily mini-291 mum, maximum, and mean temperatures; global radiation; rel-292 ative humidity; wind speed; and precipitation (the methodology293 is described in Meitner et al. 2023).294 For the Czech Republic, the selected models included295 CMCC-ESM2, EC-EARTH, GFDL-ESM4, MPI-ESM1-2-HR,296 MRI-ESM2-0, and TAIESM1. Owing to their spatial resolution,297 GCM simulations are not directly applicable for assessing local-298 ized impacts. Therefore, either systematic bias correction or299 transformation of the observed series is necessary. We applied300 the "delta method", which is a preferred approach in the Czech301 Republic for its robustness, aligning changes in the observed302 and transformed series with those in the climate model simula-303 tions. For daily resolution analysis, we employed the advanced304 delta change method (van Pelt et al.,2012;Kraaijenbrink and305 colleagues,2013), which accommodates variability changes,306 enabling differential changes in extremes relative to the mean.307 This approach also adjusts for systematic simulation errors in308 precipitation, whereas temperature transformations, which are309 linear, are unaffected by systematic error. Other meteorological310 variables (solar radiation, relative humidity, and wind speed)311 are modified by multiplication using a ratio of standard devi-312 ations for the control and scenario periods in the climate model313 simulations.314 The ADC method uses a reference period from the observed315 time series. To maintain the stationarity of the time series, the316 period 1981–2010 was selected as the reference period. On the317 basis of this reference period, four future time blocks were318 created: 2015–2044 (referred to as "2030"), 2035–2064 ("2050"),319 2055–2084("2070"), and 2070–2099 ("2085").These periods over-320 lap, allowing for smoother transitions across the projected time321 lines. 322 The final ensemble model and selected individual models 323 were applied to future climate scenario data for the Czech 324 Republic to assess the BFAI expected in the coming decades. 325 This provided additional insights into the coherence of the 326 projected results by the independent estimation approaches. 327 LOESS smoothing was used to visualize projected BFAI trends 328 across periods and the mean of climate models. 329 Results Observed trends in burnt forest areas 330 Forest fire trends were analyzed for 12 European countries using 331 BFAI data from 1990–2022 (Fig. 2). Our results present a coun332 terpoint to the broader global trends of increasing burnt forest 333 areas identified by Jones et al. (2024). Specifically, of the 12 334 European countries analyzed here, 10 countries do not show a 335 significant trend in the BFAI, while a significantly decreasing 336 trend in the BFAI is observed for Poland and Spain (Table 1;337 Fig. 3). Portugal has the highest values and variability in the 338 BFAI, followed by Spain and Mediterranean countries, includ339 ing Greece, Croatia, and Italy. As expected, the least burnt 340 areas are found in Austria and Sweden, which represent high341 elevation and Nordic conditions in Europe, respectively. 342 Model estimation of burnt forest areas 343 The machine learning analysis revealed substantial variation 344 in model performance across the 12 algorithms tested (Fig. 4). 345 The RMSE values for the calibration dataset ranged from 0.42 346 (Bagged Classification and Regression Trees) to 1.01 (Support 347 Vector Machine with linear kernel), whereas the verification 348 RMSE values ranged from 0.87 (Random Forest) to 0.94 (Bagged 349 Trees and Elastic Net). The close agreement between the calibra350 tion and verification performance across most models indicates 351 good generalization capability without substantial overfitting. 352 The performance of the tree-based methods was consistently 353 strong, with the Random Forest method achieving the best ver354 ification performance (RMSE = 0.70 for calibration and 0.87 355 for verification). Bagged Classification and Regression Trees 356 showed excellent calibration performance (RMSE = 0.42) but 357 higher verification error (RMSE = 0.94), suggesting some over358 fitting. Extreme Gradient Boosting performed well, with RMSE 359 values of 0.93 (calibration) and 0.88 (verification), whereas 360 Cubist Regression yielded RMSE values of 0.99 (calibration) and 361 0.88 (verification). 362 Linear and flexible regression methods showed moder363 ate performance, with the Generalized Linear Model achiev364 ing RMSE values of 0.99 (calibration) and 0.89 (verification), 365 respectively, while Elastic Net regression performed similarly, 366 with 0.96 (calibration) and 0.94 (verification), respectively. 367 Multivariate Adaptive Regression Splines (MARS) showed com368 5 Cienciala et al. International Journal of Wildland Fire Fig. 3. The observed trends in BFAI for 12 European countries for 1990–2022. There is no trend in ten countries, while the BFAI decreases significantly in Poland and Spain. 6 www.publish.csiro.au/xx International Journal of Wildland Fire Table 1 Analysis of trends in the BFAI in 12 European countries (Country), showing projected tendencies (Trend), coefficient of determination (R2), and the slope and p value of the linear regression between the BFAI and Year; the p values were determined via the Mann–Kendall test. Country Trend Slope R2Lin. reg. p-value Mann–Kendall p-value AUT Stable 0.000 0.007 0.647 0.588 BGR Stable -0.041 0.016 0.487 0.768 CZE Stable -0.007 0.056 0.185 0.676 ESP Decreasing -0.140 0.191 0.011 0.007 FRA Stable -0.017 0.066 0.149 0.034 GRC Stable -0.032 0.019 0.444 0.075 HRV Stable -0.008 0.001 0.899 0.566 ITA Stable -0.065 0.065 0.153 0.057 POL Decreasing -0.040 0.191 0.011 <0.001 PRT Stable 0.147 0.006 0.657 0.525 SVN Stable -0.005 0.008 0.623 0.001 SWE Stable 0.003 0.046 0.230 0.097 parable performance, with RMSE values of 0.99 (calibration)369 and 0.89 (verification).370 The ensemble model (Fig. 5), constructed through stacked371 generalization, achieved the best overall performance, with372 RMSE values of 0.81 (calibration) and 0.90 (verification), as373 shown in Fig. 6. While the ensemble did not achieve the lowest374 individual RMSE values compared with those of some individ-375 ual models (Random Forest had a RMSE of 0.87 for the veri-376 fication dataset), it provided robust and balanced performance377 across both datasets. The ensemble approach effectively com-378 bined the strengths of the different modeling approaches, with379 R2values of 0.84 and 0.80 for the calibration and verification380 datasets, respectively.381 Analysis of the ensemble composition revealed that Cubist382 Regression contributed most heavily to the final predictions,383 accounting for approximately 35% of the ensemble weight384 (Fig. 6). The Artificial Neural Network and Support Vector385 Machine (Radial) were the second and third most influential386 components, contributing approximately 25% and 15% of the387 ensemble weight, respectively. Multiple instances of Extreme388 Gradient Boosting and Elastic Net regression also contributed to389 the ensemble, reflecting the value of model diversity in improv-390 ing predictive performance.391 The observed versus predicted plots (Figures 4 and 6) demon-392 strate that the models effectively captured the substantial varia-393 tion in the BFAI across countries. Portugal consistently had the394 highest predicted and observed values (log(BFAI) >3), followed395 by Spain and Mediterranean countries. The lowest values were396 consistently recorded in Sweden and Austria ((log(BFAI) <0),397 with the models accurately capturing this geographic gradient.398 The scatter around the 1:1 line was generally modest across399 the prediction range, although some heteroscedasticity was evi400 dent, with larger residuals at higher BFAI values, particularly 401 for extreme fire years in Mediterranean countries. 402 Model projection of burnt forest areas 403 The ensemble model was applied to project the future BFAI for 404 the Czech Republic under the SSP2-4.5 climate scenario using 405 six global climate models (GCMs) for four future periods: 2030 406 (2015–2044), 2050 (2035–2064), 2070 (2055–2084), and 2085 407 (2070–2099). The projections reveal a consistent upward trajec408 tory in the burnt forest area across most climate models, despite 409 the historical success of fire management in maintaining stable 410 burn areas (Figure 7). 411 The reference period (1991–2020) had a mean observed BFAI 412 of less than 0.15 ha/kha for the country, represented by the hori413 zontal dashed line in Figure 7. This baseline reflects the current 414 fire regime under existing climate conditions and management 415 practices. All six GCMs project increases in the BFAI relative 416 to the historical baseline, with projections ranging from mod417 est increases to more substantial increases by 2085, depending 418 on the climate model. The ensemble mean, represented by the 419 LOESS-smoothed gray curve, gradually increases from under < 420 0.15 ha/kha to approximately 0.25–0.30 ha/kha by 2085, repre421 senting an approximately 1.5–2.0–fold increase over historical 422 values. 423 Substantial differences exist among GCM projections, high424 lighting climate model uncertainty in future fire risk assess425 ments. TAIESM1 consistently projects the highest fire risk 426 across all periods, reaching BFAI values of approximately 0.4 427 ha/kha by 2085. In contrast, GFDL-ESM4 and MRI-ESM2-0 428 project more conservative increases, with 2085 values remain429 7 Cienciala et al. International Journal of Wildland Fire Fig. 4. Observed versus predicted values of log-transformed burnt forest area index (log(BFAI)) across 12 models, with their performance evaluated on both the calibration (open symbols) and verification (filled symbols) datasets. Each panel represents one model, showing predictions for all 12 countries using country-specific shapes and colors. The root mean square error (RMSE) is reported separately for the calibration and verification subsets. The dashed diagonal line indicates perfect agreement. Differences in performance across models highlight variations in generalizability and predictive accuracy. 8 www.publish.csiro.au/xx International Journal of Wildland Fire Fig. 5. Predicted versus observed values of the log-transformed BFAI (log(BFAI)) for the final stacked ensemble model. Each point represents a calibration (open symbols) or verification (filled symbols) observation across the 12 countries, with distinct colors and shapes used for country differentiation. The diagonal dashed line indicates perfect agreement between the observed and predicted values. The root mean square error (RMSE) is reported separately for the calibration and verification subsets. ing closer to current levels. The CMCC-ESM2, EC-EARTH3,430 and MPI-ESM1-2-HR models show intermediate projections431 between these extremes.432 The projected increases are not linear, with most models433 showing accelerating rates of change in the latter half of the 21st434 century. The period from 2050 to 2085 shows steeper increases435 than the 2030 to 2050 period across most GCMs, suggesting that436 fire risk may increase exponentially rather than linearly with437 continued climate change.438 Discussion Trends in burnt forest areas439 Our analysis of burnt forest area, or more specifically, BFAI440 trends across 12 European countries from 1990–2022, reveals a441 contrast to global patterns of increasing wildfire activity (Jones442 et al.,2024;Zhang et al.,2024) but supports other earlier global443 observations of burnt areas (Andela et al.,2017). While 10 of the444 12 countries showed no significant BFAI trends, two countries445 (Poland and Spain) demonstrated even significantly decreasing446 trends. Additionally, although not significant for the studied447 period, most countries (9 of 12 countries) showed a negative448 slope of the regression line, indicating a declining pattern of449 the BFAI (Fig. 3, Table 1). These findings directly challenge the450 Fig. 6. Relative contribution of individual base learners to the final stacked ensemble model. Stacking coefficients were estimated using penalized regression, where the penalty parameter controls the amount of regularization applied to the model weights. A higher penalty shrinks coefficients toward zero, reducing the model complexity and mitigating overfitting. Each bar represents the stacking coefficient assigned to a model during the blending process, reflecting its influence on the ensemble predictions. Models with higher coefficients contributed more substantially to the ensemble. assumption that climate-driven increases in fire risk automati451 cally translate to proportional increases in burnt area. 452 The absence of increasing BFAI trends, despite docu453 mented increases in fire-conducive weather conditions (VPD 454 increased significantly in six countries, dead fuel moisture 455 trends increased in four countries, and the FWI increased in 456 three countries), provides compelling evidence for the effective457 ness of European fire management systems. This apparent para458 dox demonstrates that technological and institutional improve459 ments can successfully counteract climate-driven increases in 460 fire risk, at least within the European context. These find461 ings align with previous observations in Mediterranean Europe, 462 where Turco et al. (2016) reported decreasing fire trends despite 463 increasingly arid conditions. The effectiveness of European fire 464 management has been attributed to several factors, includ465 ing improved early warning systems (San-Miguel-Ayanz et al.,466 2013), enhanced suppression capabilities (Fernandes et al.,467 2016), and strategic fuel management through prescribed burn468 ing and silvicultural treatments (Fernandes and Botelho,2003;469 Moreira et al.,2011). Additionally, socioeconomic changes, 470 including rural depopulation and agricultural abandonment, 471 have altered fire regimes across Europe (Pausas and Keeley,472 2014;Moritz et al.,2014). 473 9