Cross‐country risk quantification of extreme wildfires in Mediterranean Europe

We estimate the country‐level risk of extreme wildfires defined by burned area (BA) for Mediterranean Europe and carry out a cross‐country comparison. To this end, we avail of the European Forest Fire Information System (EFFIS) geospatial data from 2006 to 2019 to perform an extreme value analysis. More specifically, we apply a point process characterization of wildfire extremes using maximum likelihood estimation. By modeling covariates, we also evaluate potential trends and correlations with commonly known factors that drive or affect wildfire occurrence, such as the Fire Weather Index as a proxy for meteorological conditions, population density, land cover type, and seasonality. We find that the highest risk of extreme wildfires is in Portugal (PT), followed by Greece (GR), Spain (ES), and Italy (IT) with a 10‐year BA return level of 50'338 ha, 33'242 ha, 25'165 ha, and 8'966 ha, respectively. Coupling our results with existing estimates of the monetary impact of large wildfires suggests expected losses of 162–439 million € (PT), 81–219 million € (ES), 41–290 million € (GR), and 18–78 million € (IT) for such 10‐year return period events.


INTRODUCTION
Wildfires affect humans, assets, and ecosystems and can lead to extensive socioeconomic and environmental impacts (Keeley et al., 2012). Within Europe, the Mediterranean region is the most fire prone with high wildfire incidence and consequences (San-Miguel-Ayanz et al., 2020). This was bleakly illustrated by the pronounced wildfire season in 2017 with blazing fires in France and roughly 140'000 hectares (ha) burnt in Portugal (NATURE, 2017), or by the 2018 fatal fires in Greece leading to more than 100 deaths and causing major damage to the ecosystems of the susceptible Natura 2000 protected areas (San-Miguel-Ayanz et al., 2018). Not fertilize soils and control plant growth) and landscape modifications (Santín & Doerr, 2016). However, although societies and ecosystems are likely to adapt to near-normal conditions, this is arguably not the case for extreme events (Bowman et al., 2017;San-Miguel-Ayanz et al., 2013;Tedim et al., 2018). Rather, evidence shows that particularly large wildfires are linked to severe disturbances and losses and are the cause of the bulk of social, economic, and adverse environmental impacts (Evin et al., 2018;Gill & Allan, 2008;Mendes et al., 2010). The purpose of this study is to characterize the spatiotemporal distribution and dynamics of extremely large wildfires in Mediterranean Europe, as well as to quantify and compare their risk probabilities across countries. Since our interest lies in the risk quantification of rare or extreme events, we model the probabilistic structure of the commonly heavy-tailed right tail of the wildfire BA density distribution (Beverly & Martell, 2005;Hernandez et al., 2015;Scotto et al., 2014) applying extreme value theory (EVT). A series of commonly known variables that potentially influence the production of large wildfires, such as the Fire Weather Index (FWI), population density, land cover type, and seasonality, are included as covariates to evaluate potential conditional probabilities. Employing the analytical tools provided by EVT enables to extrapolate wildfires of potentially unobserved size based on the European Forest Fire Information System (EFFIS) BA data set from 2006 to 2019, and thus, to quantify the risk of country-level extreme wildfires. Furthermore, we convert our estimates into rough monetary losses using figures from the existing literature to facilitate the potential application of our estimates to policy decisions.
EVT has been proven to be a suitable inferential tool for wildfire size risk quantification (Hernandez et al., 2015;Holmes et al., 2008), and has been applied globally (Jiang & Zhuang, 2011;Keyser & Westerling, 2019). For Mediterranean Europe, Evin et al. (2018) evaluate the risk of large wildfires in France conditional on a new fire policy introduced in 1994. Moreover, several studies quantify and compare regional wildfire risk and regimes in Portugal (De Zea Bermudez et al., 2009;Mendes et al., 2010;Scotto et al., 2014), which is unsurprising given the country bears the highest wildfire prevalence within Mediterranean Europe (Turco et al., 2019). However, to the best of our knowledge, ours is the first study to use EVT to perform a cross-country quantification of wildfire risk. Our contribution is threefold. First, we merge high-quality homogenized and up-to-date geospatial data sets for the European Mediterranean region. Second, we perform a country-level analysis of extremely large wildfires in Mediterranean Europe, and third, we compare the estimated risks across the region.
The remainder of this study is organized as follows. Section 2 describes the data sources followed by Section 3 outlining the methodology underpinning the extreme value analysis. Section 4 summarizes the results and derives monetary losses by matching our estimates with economic loss figures from the existing literature. The findings are subsequently discussed in Section 5, before Section 6 concludes.

Burned area
We use a high-quality BA spatial data product compiled by the Joint Research Centre (JRC) and provided by the EFFIS. 1 It is the primary source of harmonized data on wildfires in Europe, and thus enables a sound cross-country comparison. The data product is derived from semiautomatic classification of daily processing of Moderate Resolution Imaging Spectroradiometer (MODIS) satellite imagery at 250 m spatial resolution. The definite fire perimeters are refined through visual image interpretation and systematically collected fire news from various media. The data set includes fires larger than approximately 30 ha, and contains information on the initial date, country, province, place, as well as the BA polygons. 2 We model the extreme BA conditional on the covariates described hereafter. For the full list of covariates refer to Table A1 in Appendix A in the Supporting Information.

Fire Weather Index
Weather conditions are a major driver of wildfire events and are commonly applied to construct fire danger indices (Bedia et al., 2014;Krawchuk et al., 2009;Sousa et al., 2015). We employ the FWI component of the Canadian Forest Fire Weather Index System as a proxy for meteorological conditions incorporating temperature, wind speed, relative humidity, and precipitation. Providing a homogeneous numerical rating of relative fire potential resulting from the combination of the two fire behavior indices, namely, the Buildup Index and the Initial Spread Index (Van Wagner & Pickett, 1985), the FWI has become a reference index for European fire danger maps produced by the JRC (Camia et al., 2008). We use a high-resolution calculation developed by Natural Resources Canada based on the European Centre for Medium-Range Weather Forecasts ERA5-HRES 3 reanalysis product presented in McElhinny et al. (2020). To account for the effect of interseasonal drought, we use the FWI version derived from the overwintered Drought Code with a spatial resolution of 31 km (0.28 • on a reduced Gaussian grid). We spatially join the centroid of every wildfire polygon to the closest grid cell of the FWI data set and extract (i) the daily FWI values 1 month prior, as well as 1 week after the initial date of the fire, and (ii) the daily FWI values of the respective year of the fire. Using (i), we create the variables FWI on the initial date (FWI_InitDat), the mean FWI of the month prior to the initial date (FWI_MP), the mean FWI of the week prior to the initial date (FWI_WP), and the mean FWI of the month prior until the week after the initial date (FWI_MP_WA). We employ (ii) to estimate the annual mean as well as the 0.5, 0.9, 0.95, 0.99 quantiles (FWI_Mean,FWI_q0.5,FWI_q0.9,FWI_q0.95,FWI_q0.99) of the FWI for the corresponding year of the fire incidence.

Population density
Population density has gained widespread attention for its role as an ignition source, as a facilitator of suppression efforts, and as a factor that captures impact-related importance (Fernandes, 2019;González-Cabán, 2009;Lankoande & Yoder, 2006;Pechony & Shindell, 2010). To proxy population density near wildfires, we use the Oak Ridge National Laboratory's LandScan 4 annual global population distribution data provided at approximately 1 km (30'') spatial resolution. The raster data representing the ambient population distribution are based on remote sensing imagery analysis techniques, and demographic and geographic data.
We create approximately 4 km buffers 5 around the centroid of the polygons and calculate the mean population density in counts per square kilometer of the respective LandScan year denoted by the variable Pop_4km.

Land cover type
The 2006, 2012, and 2018 versions of the Copernicus' CORINE land cover 6 are employed to categorize the EFFIS perimeters of BA to evaluate a potential correlation between land cover type and the distribution of large fires. The CORINE land cover information is derived from satellite data 7 using a minimum mapping unit of 25 ha and consists of an inventory of 44 land cover classes. We extract the dominant land cover type for each EFFIS BA polygon considering the latest version of the CORINE land cover data with respect to the initial date of the observation. We further reclassify the most prevalent land cover types for each country with regard to the extreme wildfire observations. For an overview of the dominant land cover types for each country, see Table A2 of Appendix A in the Supporting Information. In the conducted analysis, types I to III are incorporated as indicator variables. A fourth indicator variable named Type_Other is created where none of the three types applies.

Point process (PP) using maximum likelihood estimation
The foundation of the PP approach regarding extremal processes was originally introduced by Pickands (1971), and applied to environmental processes by Smith (1989). The PP approach is particularly suitable as it uses data efficiently and can easily be adapted to include temporal or covariate effects (Coles, 2001). We apply a nonhomogeneous PP model to simulate the occurrence (i.e., frequency of exceedance) and intensity (i.e., excess) of a value of BA above a chosen threshold.
Let X i be a series of independent and identically distributed (i.i.d.) random variables representing wildfire BAs, and N n = {( i n+1 , X i ) : i = 1, n} be a sequence of PPs. Then, given a sufficiently large threshold u, on regions of the form [0, 1] × (u, ∞), the PP N n is approximately a Poisson process with the intensity measure Λ(A) shown in Equation (1) on a set of the form A = [t 1 , t 2 ] × (u, ∞): The interval (t 1 , t 2 ) on the abscissa is a subset of [0,1] and n y denotes the number of years of observations so that events in nonoverlapping subsets of [0, 1] × (u, ∞) are independent and the estimated parameters , , and correspond to the generalized extreme value (GEV) distribution. It envelops three types of limit distributions, which are uniquely defined by the shape parameter . The Fréchet distribution ( > 0) is characterized by a heavy tail, the Gumbel distribution ( = 0) exposes an exponential decay of the tail, whereas a Weibull limit distribution ( < 0) has an upper bound. In general, a heavier tail implies that the probability of an "unexpected" event is larger, while the location and the scale parameters relate to the mean and spread of the distribution, respectively. For greater detail on the GEV see Appendix B in the Supporting Information.
Following Coles (2001), the model parameters are estimated by maximizing the likelihood function where i is one if the realization of X i > u, and zero otherwise. The first part of the likelihood expression entails the contribution of the number of fire events (occurrence) characterized by the Poisson distribution with mean Λ{[0, 1] × [u, ∞)}. The second part shows the excess contribution of the observations (intensity), which are modeled as generalized Pareto distribution (GPD). is adjusted as * = (u) − u, so that the scale parameter is independent of the threshold. The cumulative GPD is given by Equation (3): where 1 + ( z−u * ) > 0, z − u > 0, and * > 0.
Given that X has a GPD, the distribution of the rescaled random variable z∕ * is independent of * (Katz et al., 2005).
We perform the numerical optimization using the R package extRemes (Gilleland & Katz, 2016) to estimate Equation (2).

Model assumptions
The theoretical justification for using a PP characterization of extremes is predicated on the assumptions of (i) unbiased threshold choice, (ii) stationarity, and (iii) independence of the excesses. Regarding (i), too low a threshold leads to a bias potentially violating the asymptotic basis of the model. If the threshold chosen is too high, the reduction of data points leads to high variance. We determine the individual countries' thresholds using the threshold diagnostic tools provided in the R package extRemes. They are based on the following rationale. Let the excesses over a threshold u be defined as y = x − u. Recalling from Section 3.1 that these excesses follow a GPD, this also holds true for all y > 0 of a threshold v > u with As a consequence, Equation (4) can only be satisfied if This implies that for a sufficiently high threshold, both the shape parameter and the modified scale − u are independent of the threshold and need to be stable. Besides plotting the shape and modified scale parameters individually, the mean value of the excesses y over u can be plotted against u, which is known as the Mean Residual Life (MRL) plot (Coles, 2001). The GPD is deemed to fit the data well when a straight line starting from the selected threshold can be fitted within the confidence bands of the MRL plot, and thereby indicating a stable distribution. In practice, the visual interpretation of the MRL plot, as well as the individual parameter plots, is somewhat subjective. Thus, we additionally consider the threshold selection suggestions provided by the automated Bayesian leave-oneout cross-validation approach, which compares the extreme value predictive performance resulting from each of a set of thresholds. This approach was first introduced by Northrop et al. (2017) and is implemented in the R package threshr (Northrop & Attalides, 2020). As this approach is only applicable for independent observations, we compare outcomes in an iterative process with the estimation of the extremal index as a measure of dependence. While Equation (2) implicitly assumes stationarity of the GEV parameters, we also estimate nonstationary models where , , and are conditioned on various functional forms of the covariates described in Section 2 (as well as on seasonality variables) in order to assess assumption (ii). Equation (5) serves as an example of modeling a nonconstant linear loca-tion parameter dependent on the mean FWI of the month prior to the initial date FWI_MP: The evaluation of the nonstationary models is based on the Akaike information criterion (AIC), the Bayesian information criterion (BIC), and, for nested models, on the likelihood ratio test. 8 We systematically model the location, shape, and scale parameters individually and combined starting with linear functional forms of the respective parameters. Whenever a model shows an improvement over the stationary model, we explore more complex functional forms (i.e., quadratic and interactions). In cases where the parameter confidence intervals (CIs) could not be estimated via the delta method, 500 iterative bootstraps with replacement were applied to evaluate the parameter significance.
As a means to examine the independence assumption (iii), the degree of dependence is explored using the extremal index ∈ (0, 1] suggested by Ferro and Segers (2003), which is defined as: where T i denotes the length between excesses (interexceedance time). A value of the extremal index = 1 implies complete independence, whereas → 0 indicates perfect dependence. In case the extremal index suggests a violation of the independence assumption, the data can be declustered to filter the dependent observations.

Model fit
In order to assess the fit of the selected model, we implement two common diagnostic plots incorporated in the extRemes package. First, we avail of the Z-plot following Smith and Shively (1995). Let Zk be the Poisson intensity parameter integrated from exceedance time k − 1 to exceedance time k (starting the series with k = 1). The Z-plot then determines whether this random variable Zk is independent exponentially distributed with mean one, which corresponds to the observations lying on the diagonal. Second, we plot the kernel density functions of the observed data versus the modeled distribution. For the particular case of the PP characterization of extremes, the density of the calculated data block maxima is compared to the PP model with respect to the equivalent GEV.

Return levels (quantiles)
Harnessing the estimated probabilities associated with extremes, the interest is typically focused on providing estimates of the upper quantiles of the modeled distribution functions. Specifically, the return level of an extremely large fire, defined as z p , which is associated with a return period of 1∕p, embodies a tangible outcome. It is equivalent to the (1 − p)th quantile of the corresponding modeled distribution by the PP representation of extremes. As the PP approach combines the Poisson distribution parameter with the GPD, the return level z p is obtained by setting the cumulative distribution function Equation (3) equal to the desired quantile 1 − p. Solving for z (for a probability p) leads to Equation (7) (Coles, 2001): where the return level z p denotes the BA level that is expected to be exceeded in any given year with probability p.

Economic valuation
A transformation of the informational content of the BA return-level estimates into economic values could arguably be beneficial supporting policy decisions. Our approach in this regard is to multiply associated per ha monetary losses with the expected BA, as suggested by Holmes et al. (2008). To this end, we resort to existing studies either providing explicit per ha loss estimates or we calculate per ha values by combining information on total BA with total loss estimates. To facilitate a comparison over space and time, spatial values are harmonized in hectares, and monetary values are inflation adjusted and expressed in 2020 euros using the 2020 monthly average exchange rate (US$ = 0.87€). Values in U.S. dollars are deflated based on the not seasonally adjusted urban Consumer Price Index CPI. 9 Table 1 provides a summary of the economic impact aspects that have been included in each study as well as of the inflation-adjusted €per ha monetary values for the five papers included in the table. Note that the calculated €/ha losses are highly dependent on the estimation method used, the type of damage and losses included, and on the specific situation of the fire (season) that is studied in the article. Therefore, the figures derived from the multiplication of our return levels with estimates from existing publications need to be interpreted with caution.
More specifically, two studies have a European context. The first estimates we derive are from a comprehensive report for Mediterranean forests by Merlo and Croitoru (2005), who provide figures in 2001 prices that encompass country-9 https://www.bls.gov/cpi. specific estimates of 884€/ha (GR), 1'480€/ha (IT), and 3'420€/ha (PT). This translates into inflation-adjusted monetary values expressed in 2020 euros of 1'228€/ha (GR), 2'004€/ha (IT), and 4'728€/ha (PT). We also include the economic impact estimates from a study of Galicia, Spain, by Barrio et al. (2007), which implements an ecosystem service approach based on assessing services that are affected due to wildfire existence. The reported monetary losses range from 2'249 to 3'162 €/ha in 2006 values. We apply the mean of this range, that is, 3'304€/ha in 2020 euros.
As suitable research on Southern Europe is limited, we also include three loss estimates derived from studies of U.S. wildfires. The Butry et al. (2001) case study assesses the Florida 1998 summer wildfires that burned a total of around 500'000 acres (202'343 ha). We apply the conservative lower bound total cost estimate of US$ 600 million (in 2001 values). Dividing the total cost by the total BA leads to an inflation-adjusted estimate of 3'801€/ha. A second study by Rahn et al. (2014) evaluates the 2003 wildfires in San Diego (United States) and reports a cost of US$ 6'500 per acre (US$ 2'630/ha in 2014 values), which is equivalent to 3'230€/ha in 2020 euros. Finally, we include the recent publication on California wildfires by Safford et al. (2022), who investigate the extraordinary 2020 fire season. The authors estimate losses of US$ 19 billion for a historical record of 1.74 million ha BA, which is equivalent to 8'717€/ha in 2020 euros. 10

Summary statistics
Country-level BA summary statistics of the EFFIS BA data product are presented in Table 2. The single largest fire in the data set that burnt 67'521 ha occurred in October 2017 in Portugal. Greece exhibits a comparably lower frequency of fires but has the largest mean, median, 75 percentile, and 90 percentile BA values. The highest annual wildfire incidence is recorded in Italy. Figure 1 presents the log transformed wildfire BA observations from 2006 to 2019 at the country level. The data show a slightly decreasing tendency in BA for Spain, France, Italy, and Greece, and no clear trend in Portugal. However, a slightly different picture emerges from Figure A1 of Appendix A in the Supporting Information when we focus on BA extremes, herein defined as the wildfires that exceed the selected country-level threshold. While, once again, we observe decreasing BA trends for France, Italy, and Greece, 10 Extensive research conducted by Wang et al. (2021) estimates the economic losses for the 2018 wildfire season in California that include indirect losses and suggests total wildfire damages were in the region of US$148.5 billion for a total BA of 7'700 km 2 . This leads to a per ha loss estimate of 165'467 €/ha in 2020 euros, which is around 19 times larger than the Safford et al. (2022) estimate for the 2020 season. Given this estimate is far beyond all the other estimates, we do not use it in this analysis and only present the more conservative estimates. However, Wang et al. (2021) give some indication of how far reaching the costs of extreme wildfires are when we include indirect health costs as well as costs arising outside the affected region assuming extreme fires in Mediterranean Europe are comparable to those in California.

TA B L E 1
Included economic impacts and €/ha economic loss estimates in 2020 euros   Figure 2 displays the annual number of wildfires and total BA over the study period enabling a direct country comparison. There are few observations and little variation over the years for Greece, while the opposite is the case for Italy. The year 2017 particularly stands out for Portugal with many fire records as well as large total BA. For France, 2019 accounts for more than half of all the observations in the study period.
Regarding the correlation between BA and the covariates, there is a general tendency of a positive correlation between the BA and the mean FWI of the week prior to the fire initial date (FWI_WP). No conclusive relationship is observable between the BA and population density. Correlation plots are provided in the Supporting Information in Appendix A ( Figure A2 and Figure A3 for the association of the FWI and the population density, respectively).

Threshold selection/dependence test
The MRL plots with the final threshold choice (after considering all decision-supporting tools outlined in this section) are shown in Figure 3. Complementing the MRL plots, the individual behavior of the shape parameter and the modified scale parameter − u are analyzed but not shown. The Bayesian leave-one-out cross-validation plots are presented in Figure 4. These show a single run output and vary across different executions. The best threshold evaluated by this approach, denoted as u b , is provided below the plot whenever it proved stable over 10 consecutive runs. Table 3 provides a summary of the final threshold choices with the corresponding extremal indices , and the number of observations above the selected threshold. The excesses of Spain and Greece indicate perfect independence, while Portugal and Italy show a very high value. The lowest extremal index value is found for France, which did not improve after

Country
Threshold u Extremal index n > u (% of total)  (6) declustering. 11 The number of observations above the respective country-specific thresholds ranges from 42 to 62 and corresponds to 1.4%-6.3% of the total country-level data.

Nonstationarity
Since the stationary models are embedded in potential nonstationary models, the results of the latter are reported first. Table 4 lists all the models with an improvement of the BIC > 10 over the stationary model following Neath and Cavanaugh (2012), suggesting this threshold as "very strong" 11 The specific case of modeling the extremes with the data available for France is addressed in Section 4.3 and Section 4.4.
evidence to favor the model with the lower BIC over the competing model. For Portugal, letting the location parameter depend on the mean FWI for the month prior to the initial date of the fire (FWI_MP) leads to the best model fit. The evaluation of conditional effects for the historical excesses in Spain and Italy shows that modeling the location and the shape parameter dependent on the FWI on the reported initial date (FWI_InitDat) improves the model fit the most. None of the nonstationary models leads to any improvements of the model fit for Greece. For France, land cover type is found to be most influential in modeling the observed data. More specifically, modeling the location parameter conditional on land cover Type_I (Sclerophyllous vegetation) in a linear functional form not only proves to capture the empirical data best, but also leads to a significant positive shift of the distributional mean. However, even though modeling nonstationarity leads to an increased model fit in specific cases, with the exception of France, results do not indicate a significant modification of the modeled GEV parameters. Consequently, based on the covariates considered, the assumption of stationarity holds true in the data sample for all countries except for France. Therefore, reported probabilities of the stationary model are valid and comparable for Portugal, Spain, Italy, and Greece.

TA B L E 4
Nonstationary models with a BIC decrease > 10 sorted by decreasing AIC

Greece No model improvements
Note: (i) Whenever the functional form is indicated as "quadratic," the linear term is included as well. (ii) In cases where the parameter 95% CIs could not be estimated via the delta method, 500 iterative bootstraps with replacement were applied to evaluate the parameter significance (indicated with *).

Model selection/model fit
In light of no significant parameter changes modeling the extremes conditional on the implemented covariates for Portugal, Spain, Italy, and Greece, we base the subsequent model evaluations and estimations on the stationary model. Not only does the distribution of the historical extremes for France show dependence and therefore violate the stationarity assumption, but the extremal index in Table 2 also indicates higher dependence of the excesses than is the case for the other countries. Thus, we exclude France from subsequent analysis. Evaluating the model fit using Z-plots depicted in Figure 5, it is evident that all observations lie well within the 95% confidence bands and that there is arguably a good model fit for Portugal, Spain, and Italy, and a moderately good fit for Greece. A similar conclusion can be drawn from Figure 6 plotting the kernel density functions of the empirical against the modeled data. Once again, the observed data are very well modeled for Italy, and fairly well for Portugal and Spain. For Greece, the modeled data, in contrast, captures the empirical data relatively less well.

Parameter estimates
The three GEV distribution parameter estimates by country are shown in Table 5. The largest point estimatêis found for Portugal, followed by Spain, Greece, and Italy, respectively. Although the center of the distribution is larger for Spain than for Greece,̂indicates that the spread of the distribution is wider for Greece than for Spain. In general, we observe extremely small CIs for Portugal for the location and scale parameters. The largest shape parameter valuê, and thus the heaviest tail, is estimated for Greece followed by Portugal and is larger than 0.5 indicating that although the mean is finite, the variance is infinite (Katz et al., 2005). 12 The point estimates of the shape parameter for Spain and Italy are fairly similar. However,̂is insignificant for Italy. On that account, the main difference in the distributions of the extremes comparing the individual countries is that Portugal, Greece, and Spain have a significantly positive shape parameter indicating a Fréchettype limit distribution, while the excesses for Italy follow the Gumbel-type limit distribution.  Figure 7 show the distribution of the observations within the tail. Essentially, the limit distributions found for all the Mediterranean countries have no upper bound (i.e., the extremes are not converging to a specific value). Furthermore, the return level plots enable a better understanding of the different limit type distributions in a graphical fashion. In particular, the distinction between the Gumbel-type distribution found for the extremes in Italy versus the Fréchet-type distributions for the other countries is distinctly visible. As the x-axis is log-transformed, the return-level plot reflects Gumbel-type distributions characterized by an exponential decay of the tail as a straight line, while the Fréchet-type distributions manifest as convex shapes.

Return levels and probabilities of exceedance
We find that all the observed events lie within the bootstrapped CIs. Furthermore, the smallest confidence bands at the 95% level are observed for Spain indicating high certainty of the point estimates. In contrast, the largest CIs are apparent for Italy suggesting a wide range of potential outcomes within the 95% CI.
Looking at the largest wildfire for each country (max value in Table 2), an event of such size or larger is expected to occur, on average, once every 16 years for Portugal with an annual occurrence probability of 1.9%. The calculated yearly probability for the largest observed fire in Spain is 1% and has a return period of about 18 years. In Italy, the maximum BA value is expected to be exceeded once in every 23 years with an annual probability of 2.6%, and the largest BA value for Greece, is estimated as an approximately 16-year event F I G U R E 7 Return-level plots with bootstrapped CIs on the 95% level with a yearly probability of occurrence of 1.8%. Figure 8 overlays the individual country return-level plots to facilitate a cross-country comparison of the extremal BA distribution. We find the highest risk for extremely large fires for any given return period for Portugal and the lowest for Italy. Comparing Greece and Spain, a higher risk for large BAs emerges for Spain for low return periods (approximately < 3 years) but above this threshold, the return levels are distinctively larger for Greece.
A similar picture emerges when overlying the individual country-level BA thresholds that are exceeded in any given year with corresponding probabilities shown in Figure 9. The annual probability for extremely large fires decreases fastest for Italy and slowest for Portugal. The rate of the yearly probability decrease is comparably close for Spain and Greece with the parameter point estimates only differing marginally.

Economic valuation
Combining our results with the economic loss figures in €/ha leads to expected return period-specific economic losses pre-TA B L E 7 Range of country-level economic loss estimates for specific return levels (rl) in million € (in 2020€) sented in Table 7. Allowing a comparison of the individual publications' loss calculations, Figure 10 graphically displays the economic loss estimates for wildfires that are expected to occur, on average once every 20 years. Recall from Table 1, that while the estimates by Butry et al. (2001), Rahn et al. (2014), and Barrio et al. (2007) are relatively close, the country-specific €/ha estimate based on the figures in Merlo and Croitoru (2005) is lower than the other three for Italy and Greece, and higher for Portugal. The latest study conducted by Safford et al. (2022) clearly stands out with a distinctively larger loss estimate value.
In addition to providing economic loss estimates for specific return periods resulting from the extreme value mod-F I G U R E 8 All country return-level plot F I G U R E 9 Country-level BA exceedance probabilities in any given year eling, we also show cost estimates based on the single largest observed wildfires in the study period for each country. Hence, the cost estimates come from multiplying the maximum events in Table 2 by the corresponding €/ha estimates derived from the existing literature. The largest wildfire event leads to an economic loss estimate of 218-589 million € for Portugal, 105-283 million € for Spain, 23-101 million € for Italy, and 56-399 million € for Greece for a specific event of that magnitude.

Implications
With the quantification of country-level risk of extreme wildfires, we are able to contribute to the empirical evidence to information-based decision making regarding forest management for various stakeholders. Providing reliable estimates of return periods arguably has important implications for government agencies looking to adjust budget planning for fire prevention measures and suppression spending. Furthermore, the quantification of large fire risk through return levels can provide useful information for landowners regarding long-term investment and forest management choices, or for other institutions such as reinsurance companies. Moreover, the knowledge of the wildfire risk could also be used to increase awareness and thus may affect decision making at the individual level (i.e., location choices, property protection measures, investment in insurance associated with wildfire damage). Converting the return-level estimates of extreme wildfires in Mediterranean Europe to monetary values, as we did here, arguably provides an important tool for policy-related cost-benefit analyses. For example, the associated monetary values with a return period event can assist a government in the budget allocation of both fire prevention and suppression spending by comparing their expenditures with the expected losses F I G U R E 1 0 Country-level economic loss estimates for the 20-year return period particularly for extremely large wildfires over a specific time period.
Examining the specific results, it is insightful to first reflect on the implications of the different distributions of extreme wildfires estimated for the individual countries in our analysis. Most importantly, we find that these rare events follow a Fréchet-type distribution for Portugal, Spain, and Greece. This is in line with regional estimates within Portugal by De Zea Bermudez et al. (2009) andScotto et al. (2014). Out of the three limit-type distributions, the Fréchet distribution has the heaviest tail indicating that the probability of rare events is much higher than commonly perceived (i.e., "extreme" wildfires are not as surprising). However, although the point estimate of the shape parameter is very similar for Spain and Italy, it is not found to be significantly positive for Italy, implying that the respective extremes follow a Gumbel distribution characterized by a lighter tail than for the other Mediterranean countries. Thus, extremely large wildfires are expected to occur less often in Italy than in Spain. Overall, we find the largest point estimate of the shape parameter of 0.58 for Greece, followed by a value of 0.52 for Portugal. This indicates that the probability of extremely large wildfires is highest in Greece when only the shape parameter of the extreme event distribution is considered (i.e., excluding the mean and the spread of the distribution). Notably, both the Fréchet and the Gumbel distribution do not converge to an upper limit but are unbounded. As a matter of fact, the extremes, and thus the associated losses, characterized by the Fréchet distribution, are limitless.
The return-level results derived from the inclusion of all the three parameter estimates (location, scale, and shape) indicate the highest risk of extremely large wildfires across all evaluated return periods in Portugal and the lowest risk in Italy. Comparing Greece and Portugal, the return level for up to about 3-year events is higher in Spain, but for any return period above, it is found to be higher for Greece. For instance, the individual country return levels for 10-year return period events are 50'338 ha ( Our data do not suggest that the FWI, which captures relative fire danger, affects the distribution of wildfire occurrence or their magnitude. Ideally, we would have a much longer time series, which would make it possible to detect climatological changes. This means our results should be interpreted with care. In this regard, it must be pointed out that while many climate change projections suggest that Southern Europe faces an increasing risk in extreme wildfires (Bowman et al., 2017;De Rigo et al., 2017;Turco et al., 2018), Batllori et al. (2013) indicate that fire activity predictions can be highly divergent particularly regarding precipitation-related variables. Notwithstanding the wildfire risk driven by future climate conditions, evidence also suggests that the risk associated with human exposure may increase especially with projected population growth in fireprone regions (Knorr et al., 2016;Turco et al., 2019). Even though we did not find evidence supporting a time trend in our study, it is crucial to continue efforts to better understand the risk associated with wildfires for Mediterranean Europe. Going forward, more comprehensive and harmonized data are needed to evaluate future extreme wildfire risk scenarios incorporating climatic and demographic components as well as more detailed information on the individual fires (e.g., duration, severity, ignition point, cause) in order to distinguish which factors have the potential to influence the extremely large fires.

Limitations
Although predictions of events not actually observed in the historical data are common with the use of EVT methods, we need to emphasize that our estimates are based on data for a particularly short time period of 14 years. In this regard, encounter probability 13 suggests that the probability 13 The encounter probability P e = 1 − (1 − 1 T ) n is the likelihood of observing a T-return period event within a specific time period denoted by n.
of observing a 5-year event given our data is approximately 96%, a 10-year event 77%, a 20-year event 51%, and a 50-year event 25%, and only our 10-and 20-year event estimates are based on events witnessed in the sample period. The short time period of data may also play a role in the nonstationarity results. Although modeling the threshold excesses conditional on factors potentially influencing the distribution of extreme wildfires does lead to an improved capture of the empirical data for all countries but Greece, none of the models significantly changed the extremal distribution. Whereas it is not given that the included variables would lead to a change in the distribution of extremely large wildfires with a prolonged time period, a potential underlying dependence on factors driving or affecting extreme BA is more difficult to detect in shorter periods of analysis.
The coupling of our estimated BA return levels with existing economic loss figures also comes with strong caveats, particularly with regard to regional and temporal transfers of monetary estimates, as well as through the distinct study designs incorporating disparate economic variables in the respective loss calculations. For example, only Barrio et al. (2007) and Butry et al. (2001) include any estimation of wildfire-related health costs, which are of significant magnitude and thus, of rising concern as pointed out in Black et al. (2017). Furthermore, even though Merlo and Croitoru (2005) address country-level estimates of indirect use, option, bequest, and existence values of forests in general, they are not applied to the BA scenario and only the estimate provided by Safford et al. (2022) includes ecological (vegetation and wildlife) damage. Moreover, as the monetary valuation of indirect costs poses great challenges, besides the impediment imposed by oftentimes limited data availability particularly driven by methodological restraints, many loss calculations focus on direct impacts. However, the indirect costs are likely to exceed the reported costs as argued in CCST (2020). Thus, our calculations are conservative and the considered losses are likely to represent only some fraction of the actual economic impact. Nevertheless, our results still provide some indication of the serious implications wildfires have for many other sectors that they can reach far beyond the commonly assessed impacts.
In terms of examining the role of covariates in potential changes in the distribution of extreme wildfires, the FWI may not be the most suitable variable to capture those. Jiménez-Ruano et al. (2019) conclude that although the FWI provides useful information regarding seasonal variability and nearfuture trends, it is not necessarily the most advisable index to detect long-term trends. In this regard however, Pérez-Sánchez et al. (2017) do identify the FWI as the most suitable index for fire-risk ignition and spreading in semiarid areas such as the Iberian Peninsula. Likewise, De Rigo et al. (2017) point out that the FWI is well suited as a harmonized index over different regions for weather-driven fire danger, and Fernandes et al. (2016) observe that particularly large fires exhibit stronger responses to the severity of the fire weather.
With respect to the population density covariate, Bowman et al. (2017) demonstrate that large destructive wildfires are most likely to occur in a two-sided bounded area that excludes either very sparsely or densely populated areas, and thus highlights the underlying complexity of this interdependence. On this account, our study design calculating the mean population density around the center of a BA, might arguably be an unsatisfactory way to capture this intricate relationship, given there exists one for the extremely large fires in the studied geographical area. Additionally, we do see that many of the population density values for the extreme wildfires in our data set lie within a narrow range leading to small explanatory power of the variable. Having information on the exact ignition location of a fire might contribute to the evaluation of the potential association between population density and extreme wildfire occurrence.
Regarding the land cover type, we face a slightly different problem as it is implemented as a categorical variable. As we look at the extremes, we focus on solely on 42-62 observations for each country, and thus, we need to strictly limit the number of categories to forego having only very few wildfires for each of those. Therefore, we categorize the extreme wildfires into four main country-specific land covertype classes, and thereby sacrifice some of the specificity. Comparable drawbacks arise from the categorical covariates capturing seasonality as the extremely large wildfires are assigned to one of the four seasons. However, in this case, it is less a problem of simplification but rather one of unequally distributed observations per category, particularly as "off-wildfire season" categories arguably contain very few observations.
For the specific case of France, we note that compared to the other countries the data are more challenging to work with. Although it is typical for all the Mediterranean countries that certain years stand out with more severe fire seasons, this is particularly pronounced for the wildfire records in France with more than half of the observations coming from 2019. This in turn leads to a comparably high dependence of the extreme observations as many of the largest wildfires are recorded in a single year. Furthermore, in contrast to the other European Mediterranean countries, France geographically expands much further North, and is thus characterized by more diverse land cover types. Hence, the finding that a specific land cover type, namely, Sclerophyllous vegetation, leads to a positive shift in the distribution of the extremes might indicate that this is the vegetation type most dominant at the Mediterranean coastline and may be correlated with extreme wildfires. 14

CONCLUSION
In this article, we assemble a high-quality homogeneous upto-date geospatial data set for Mediterranean Europe and perform a cross-country risk analysis of extreme wildfires defined by BA. Although modeling a variety of covariates with the potential to affect the extremal distributions, we find no evidence for nonstationarity in the observed study period. Furthermore, the threshold excesses for France in the data set do not fulfill underlying assumptions to carry out a sound EVT analysis, and are thus only included in the descriptive part. In our results, we find the highest risk for extremely large wildfires in Portugal, followed by Greece, Spain, and Italy. We estimate the return levels for 5-, 10-, 20-, and 50-year return period events and combine our outcomes with the existing literature on economic costs. The robust estimation of extreme wildfire events underlying an evidencebased risk assessment is arguably beneficial for governmental bodies, reinsurance institutions, landowners, and residents in wildfire-prone areas providing support in information-based decision-making processes. We emphasize the need to build international homogeneous comprehensive databases with high spatial and temporal resolution regarding wildfire occurrence (ideally including point of origin, duration, and cause) but also dedicated to associated measures such as prevention and suppression spending, as well as individual fire event impact on ecosystems, infrastructures, properties, and people. Accompanying the extensive WUI with exposed communities particularly at the highly populated coastal areas of Southern Europe and vulnerable ecosystems across the region, extreme wildfire events continue to pose a substantial environmental hazard for Mediterranean Europe in the future.

A C K N O W L E D G M E N T
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement Innovative Training Networks (ITN) No. 860787.