Wildland fire potential outlooks for Portugal using meteorological indices of fire danger

Portugal is recurrently affected by large wildfire events that have serious impacts at the socio-economic and environmental levels and dramatic consequences associated with the loss of lives and the destruction of the landscape. Accordingly, seasonal forecasts are required to assist fire managers, thus contributing to alter the historically based purely reactive response. In this context, we present and discuss a statistical model to estimate the probability that the total burned area during summer will exceed a given threshold. The statistical model uses meteorological information that rates the accumulation of thermal and vegetation stress. Outlooks for the 39-year study period (1980–2018) show that, when the statistical model is applied from 26 May to 30 June, out of the six severe years, only one year is not anticipated as potentially severe and, out of the six weak years, only one is not anticipated as potentially weak. The availability of outlooks of wildfire potential with an anticipation of up to 1 month before the starting of the fire season, such as the one proposed here, may serve to provide clear directions for the fire community when planning prevention and combating fire events.

1 Introduction

Portugal is regularly affected by very large and destructive wildfires that represent a serious threat to human lives and to the territory, and have very strong and adverse impacts at the social, economic, ecologic and environmental levels that include human casualties, the destruction of homes and other structures, damages to forests, agricultural fields, shrublands, and livestock, changes in the landscape and emission of greenhouse gases (Costa et al., 2011; Pereira et al., 2011). According to the European Commission technical report about forest fires in Europe (San-Miguel-Ayanz et al., 2018), from 1980 to 2017, Portugal accounted for 26 % of the total burned area in the five southern member states (Portugal, Spain, Italy, Greece and France). However, this figure increases to 42 % when restricting to the last decade. The extreme year of 2017 deserves a special emphasis with Portugal accounting for 59 % of the burned area in the five southern member states, despite representing 6 % of the total area of the countries.

As in Mediterranean Europe, fire activity in Portugal involves complex interactions among climate, vegetation and humans (Lavorel et al., 2006; Costa et al., 2010). Persistent warm and dry conditions followed by heat spells in summer provide the optimal meteorological background to the onset and spread of large fire events (Pereira et al., 2005; Trigo et al., 2006; DaCamara et al., 2014) that take place in a landscape of cumulated biomass and increased fuel connectivity as a result of agricultural abandonment and forest expansion (Pausas and Fernández-Muñoz, 2012; Fernandes et al., 2014; Viedma et al., 2015; Oliveira et al., 2017). In addition, Portugal has the highest density of ignitions in southern Europe, with human presence and activity being the key drivers of ignitions, most of them associated with land management practices and inadequate use of fire (Catry et al., 2009). Finally, the impacts of climate change cannot be disregarded since a future warmer climate will instigate a larger number of severe wildfire episodes (Flannigan et al., 2013; Sousa et al., 2015). In this context, the catastrophic fires of June 2017 in Portugal (with 65 fatalities) associated with the strongest heatwave ever observed over Iberia in June provide a stark reminder of the growing likelihood of these events under the current warming climate (Sánchez-Benítez et al., 2018).

The magnitude of the problems related to fire activity in Portugal has motivated the scientific and technical communities to cooperate with forest managers, firefighters and the fire community with the aim of improving the understanding of the fire regime in Portugal and develop better tools that will help prevent and mitigate the impacts of severe fire events (Collins et al., 2013; Oliveira et al., 2016).

Outlooks of wildfire potential for different regions of the globe are usually based on statistical approaches (e.g. Gudmundsson et al., 2014; Turco et al., 2017) or, alternatively, using fire weather predictions a few months in advance based on dynamical seasonal forecasts by atmospheric circulation models (Anderson et al., 2007; Turco et al., 2018). For instance, the European Forest Fire Information System (EFFIS), one of the components of the Emergency Management Services in the EU Copernicus programme, is currently disseminating maps (as experimental products) of long-term seasonal forecasts of temperature and rainfall anomalies based on the European Centre for Medium-Range Weather Forecasts (ECMWF) Seasonal Forecasting System (System 4). However, as pointed out by Bedia et al. (2018), although there is significant skill in predicting 1 month ahead for above-average summer fire weather in some parts of south-eastern Europe, skill is in general quite poor elsewhere.

As an alternative, the aim of this paper is to describe and discuss a model that allows outlooks of wildfire potential in Portugal during the fire season up to 1 month in advance based on a statistical approach that integrates information about meteorological fire danger before and during the fire season. The rationale is that soil moisture deficit and drought have an impact on the increased frequency and can amplify the magnitude of hot summer extreme events in the Mediterranean (Vautard et al., 2007; Hirschi et al., 2011). Thus, persistent warm and dry conditions along the pre-fire season induce thermal and water stress on vegetation, making the landscape more prone to the occurrence of very severe fire episodes. At the same time, these conditions increase the likelihood of heatwave spells that steer the onset and propagation of large fires (Gudmundsson et al., 2014; Turco et al., 2017). The procedure involves three steps. First we set up a null model of burned area (BA) during the fire season (where no meteorological information is incorporated); then we set up a diagnostic model of BA that incorporates meteorological information during the pre-fire and the fire seasons; finally we set up a prognostic model of BA that only incorporates meteorological information before the fire season. We then show that both the diagnostic and the prognostic models perform better than the null model and that the loss in performance of the prognostic model (when compared to the diagnostic one) is relatively small, namely in respect to its capacity to anticipate fire seasons characterized by high amounts of BA. Section 2 provides a description of data and methods, and results are presented in Sect. 3. Finally, discussion of results and conclusions are drawn in Sect. 4.

2 Data and methods

The area of interest is defined as the territory of Portugal and the study covers the 39-year period from 1980 to 2018. Data of BA consist of yearly amounts of cumulated BA in July and August, hereby referred to as the fire season. BA data are derived from the official Portuguese Rural Fire Database provided by the national authority for forests (ICNF). The database contains more than half a million records of fire events, with information of total burned area, date and time of ignition and extinction, and spatial location of the starting point. Details about the database are provided in Pereira et al. (2011). Cumulated BA in the fire season accounts for more than 70 % of the total burned area in Portugal (Pereira et al., 2013), and more than 80 % of extreme fire days (defined as the top 5 % in terms of radiative energy released by wildfires) occur in July and August (DaCamara and Trigo, 2018).

Information about meteorological fire danger consists of daily values of the Daily Severity Index (DSR) covering the months of April to August. DSR is an extension of the Canadian Forest Fire Weather Index System (CFFWIS) and rates the difficulty of controlling fires (van Wagner, 1987). This index has been successfully used to model BA variability in Portugal at daily and monthly scales (Calado et al., 2008; Pereira et al., 2013). DSR results from a direct transformation of the Fire Weather Index (FWI), the last of the six components of CFFWIS, according to the relation DSR=0.0272 (FWI)^1.77. This transformation weights FWI sharply as it increases so that DSR becomes more suitable than FWI to be cumulated or averaged. For each day, the six components are computed based on consecutive daily observations of meteorological parameters of the previous days, namely temperature, relative humidity, wind speed and 24 h cumulated precipitation (Wang et al., 2015). In this study, meteorological parameters consist of gridded daily values at 12:00 UTC of 2 m temperature, relative humidity, wind speed and 24 h cumulated precipitation that were obtained from the ERA-Interim reanalysis dataset (Dee et al., 2011) issued by ECMWF. As described in Pinto et al. (2018), the original ERA-Interim data were re-projected onto the normalized geostationary projection (NGP) of Meteosat Second Generation (MSG) (EUMETSAT, 1999), with an average pixel size of about 4 km × 4 km over Portugal. Daily values of DSR for Portugal were then obtained by averaging over all grid points located within the study area.

As discussed in Pereira et al. (2005), the interannual variability of BA in Portugal is modulated by two kinds of meteorological factors, namely the temperature and precipitation regimes during the pre-fire season and the occurrence of hot and dry spells during the fire season. Following Nunes et al. (2014), the period from April to August is accordingly divided into two subperiods: (1) the pre-fire season that runs from 1 April (day 1) up to 30 June (day 91) and (2) the fire season that runs from 1 July (day 1) up to 31 August (day 62).

For each year of the period 1980–2018, meteorological fire danger at day d of the pre-fire season (pfs) of the considered year is rated by index D_pfs(d), defined as the cumulative value of daily DSR since 1 April:

$$\begin{array}{}\text{(1)}& D_{\mathrm{pfs}}\left(d\right)=\sum _{i=\mathrm{1}}^d{\mathrm{DSR}}_i,d=\mathrm{1},\mathrm{\dots },\mathrm{91},\end{array}$$

where DSR_i is the value of DSR at day i. As d runs along the pre-fire season of the year considered, progressively more information is integrated in D_pfs(d) about past daily meteorological fire danger.

In turn, meteorological fire danger of the fire season (fs) of each year of the period 1980–2018 is rated by index D_fs, defined as the square root of the mean squared anomalies performed over days characterized by a positive anomaly of DSR:

$$\begin{array}{}\text{(2)}& D_{\mathrm{fs}}=\sqrt{{\displaystyle \frac{\sum _{j=\mathrm{1}}^{\mathrm{62}}H\left[A_j\right]{\left(A_j\right)}^{\mathrm{2}}}{\sum _{j=\mathrm{1}}^{\mathrm{62}}H\left[A_j\right]}}},\end{array}$$

where H[x] is the Heaviside step function ($H\left[x\right]=\mathrm{1},x>\mathrm{0};H\left[x\right]=\mathrm{0},x\le \mathrm{0})$ and

$$\begin{array}{}\text{(3)}& A_j={\mathrm{DSR}}_j-\stackrel{\mathrm{‾}}{{\mathrm{DSR}}_j}\end{array}$$

is the anomaly of DSR at day j of the considered year, which is defined as the departure of DSR from the climatological mean $\stackrel{\mathrm{‾}}{{\mathrm{DSR}}_j}$ for day j (as obtained by averaging DSR for that day over the 39-year period 1980–2018).

The temperature and precipitation regimes during the pre-fire season are quantified by D_pfs(d), where, for a given day d of the pre-fire season, large values indicate persistent warm and/or dry conditions up to day d, inducing thermal and water stress on vegetation and a shortage of water in the soil. The impact of hot and dry spells during the fire season is in turn quantified by D_fs, which is very sensitive to the occurrence, during the considered fire season, of very large positive anomalies of DSR that are usually associated with heatwaves (Gudmundsson et al., 2014; Turco et al., 2017).

For each day d, the 39 values of D_pfs(d) for the period 1980–2018 were normalized by subtracting the respective sample mean and then dividing by the sample standard deviation. Values of D_fs were also normalized by subtracting the sample mean and dividing by the sample standard deviation. Obtained normalized indices of D_pfs(d) and D_fs will be hereby denoted as ψ(d) and χ, respectively.

The following types of models are considered in this study:

1. a null model $X\sim N(X;{\mathit{\mu }}_X{\mathit{\sigma }}_X)$, where X is the considered variable and μ_X and σ_X are the mean and the standard deviation to be estimated, respectively;

2. a nested model with two covariates $X\sim N(X;a\times A+b\times B+c,{\mathit{\sigma }}_X)$, where the mean of the normal distribution linearly depends on covariates A and B, with a, b and c being parameters to be estimated;

3. the previous nested model then simplified to a nested model with just one covariate $X\sim N(X;p\times A+q,{\mathit{\sigma }}_X)$, where the mean of the normal distribution linearly depends on covariate A, and p and q are parameters to be estimated.

Using the maximum likelihood method (Wilks, 2011), estimates of parameters μ_X, σ_X, p, q, a, b and c are obtained as follows from sample $X_j(j=\mathrm{1},\mathrm{\dots },n)$ where n is the size of the sample.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\widehat{\mathit{\mu }}}_X={\displaystyle \frac{\mathrm{1}}{n}}\mathrm{\Sigma }\left(X_j\right)\text{(5)}& {\displaystyle }{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{n}}\mathrm{\Sigma }{\left(X_j-{\widehat{\mathit{\mu }}}_X\right)}^{\mathrm{2}}\text{(6)}& {\displaystyle }a={\displaystyle \frac{\mathrm{var}\left(B\right)\mathrm{cov}(X,A)-\mathrm{cov}(A,B)\mathrm{cov}(X,B)}{\mathrm{var}\left(A\right)\mathrm{var}\left(B\right)-[\mathrm{cov}\left(A,B\right){]}^{\mathrm{2}}}}\text{(7)}& {\displaystyle }b={\displaystyle \frac{\mathrm{var}\left(A\right)\mathrm{cov}\left(X,B\right)-\mathrm{cov}\left(A,B\right)\mathrm{cov}(X,A)}{\mathrm{var}\left(A\right)\mathrm{var}\left(B\right)-[\mathrm{cov}\left(A,B\right){]}^{\mathrm{2}}}}\text{(8)}& {\displaystyle }c=\stackrel{\mathrm{‾}}{X}-a\stackrel{\mathrm{‾}}{A}-b\stackrel{\mathrm{‾}}{B}\text{(9)}& {\displaystyle }p={\displaystyle \frac{\mathrm{cov}\left(x,A\right)}{\mathrm{var}\left(A\right)}}\text{(10)}& {\displaystyle }q=\stackrel{\mathrm{‾}}{X}-p\stackrel{\mathrm{‾}}{A}\end{array}$$

Here $\stackrel{\mathrm{‾}}{\left(\right)}$, var( ) and cov( ) respectively denote the mean, the variance and the covariance operators.

In the case of the null model, the Lilliefors test (Conover, 1980) is used to test the null hypothesis that the data come from a normal distribution. In the case of the nested models with covariates, the likelihood ratio test (Wilks, 1938) is used to decide on the null hypothesis that the null model is to be retained (against the alternate model under consideration).

3 Results

3.1 Null model

The time series of cumulated BA in the fire season for the period 1980–2018 (Fig. 1) presents very large interannual variability, the extremely high amounts of 2003 and 2005 contrasting with the extremely low values that were observed in 1983, 1988, 1997 and 2008. It is worth noting that 2017, the year with the largest record in total BA (circa 450 000 ha), only ranks fourth when restricting to July and August because the largest fire events took place outside of the fire season, in June and October (Sánchez-Benítez et al., 2018).

Figure 1 Time series of yearly cumulated BA (ha) in the fire season (July and August).

The decimal logarithm of cumulated BA in the fire season follows a normal distribution model (Fig. 2); the null hypothesis that the sample of log ₁₀BA comes from a normal distribution is not rejected at the 5 % significance level by the Lilliefors test (p value of 0.21). We therefore have the following null model:

$$\begin{array}{}\text{(11)}& {\log }_{\mathrm{10}}\mathrm{BA}\sim N({\log }_{\mathrm{10}}\mathrm{BA};{\mathit{\mu }}_{\mathrm{BA}},{\mathit{\sigma }}_{\mathrm{BA}}),\end{array}$$

with μ_BA=4.69 and σ_BA=0.46 as using the maximum likelihood method.

Figure 2 Normal probability plot comparing the sample of log ₁₀BA to the normal distribution. The groups of severe and weak years are marked in red and green respectively and the years are identified by the two digits next to the symbols.

Values of log ₁₀BA above percentile 80 of the model (highlighted in red in Fig. 2) and below percentile 20 (highlighted in green) present larger departures from the fitted normal distribution than the remaining years; these two groups are classified respectively as severe years (1991, 1995, 1998, 2003, 2005 and 2017) and as weak years (1983, 1988, 1997, 2007, 2008 and 2014). The remaining years (marked in black in Fig. 2) are classified as moderate.

3.2 Model with two covariates

Departures of severe and weak years from the fitted normal distribution suggest that other factors, namely meteorological ones, are playing a role in the interannual variability of BA (Pereira et al., 2005). Following the approach proposed by Nunes et al. (2014), for a given fixed day d, where d is chosen between 26 May (d=56) and 30 June (d=91), we tested an alternate model that incorporates information about meteorological fire danger in the pre-fire and the fire seasons. The starting date of 26 May was set a posteriori, as the day when the alternate model becomes statistically significant and therefore the information provided is relevant for the users. Accordingly, using Eqs. (1)–(3), we fitted to the sample of log ₁₀BA a normal distribution with the mean linearly depending on covariate ψ(d) for the considered fixed day d, and on covariate χ; i.e. we tested the following model:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\log }_{\mathrm{10}}\mathrm{BA}\sim N({\log }_{\mathrm{10}}\mathrm{BA};\\ & a\left(d\right)\times \mathit{\psi }\left(d\right)+b\left(d\right)\times \mathit{\chi }+c,{\mathit{\sigma }}_{\mathrm{BA}})\end{array}.\end{array}$$

It is worth noting that index d is fixed in Eq. (12) and is only used to identify the fixed day of the pre-fire season when covariate ψ(d) is computed. Maximum likelihood estimates of coefficients a(d) and b(d) for the considered fixed day d are then obtained according to Eqs. (6) and (7). However, parameter c in Eq. (12) does not depend on chosen day d since, according to Eq. (8), c=μ_BA for all days because ψ(d) and χ have zero mean (since they are normalized).

An alternate model with covariates ψ(d) and χ is accordingly fitted for each fixed day d, and the null hypothesis that the null model is to be retained (against the alternate nested model) is assessed using the likelihood ratio test. For all days between 26 May (d=56) and 30 June (d=91), the null hypothesis is rejected at the 5 % significance level, with the p values steeply decreasing from $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{4}}$. In turn, when d progresses along the pre-fire season, obtained values of a(d) increase whereas corresponding values of b(d) slightly decrease (Fig. 3).

Figure 3 Temporal evolution of coefficients a and b of the normal model with two covariates.

For each day d, performance of the fitted alternate model may be assessed by representing each year in space (ψχ) framed by covariates ψ(d) and χ over a background of probability of exceedance of a given fixed threshold, e.g. of the mean value μ_BA of the null model, i.e. $P_{\mathrm{exc}}\left(\mathit{\psi },\mathit{\chi }\right)=N[{\log }_{\mathrm{10}}\mathrm{BA}>{\mathit{\mu }}_{\mathrm{BA}};a\left(d\right)\times \mathit{\psi }\left(d\right)+b\left(d\right)\times \mathit{\chi },{\mathit{\sigma }}_{\mathrm{BA}}]$. Results of models fitted on 26 May (d=56), 15 June (d=76) and 30 June (d=91) are shown in Fig. 4. Severe (weak) years tend to spread over the upper right (lower left) quadrants of the space indicating that they are associated with high (low) values of both ψ(d) and χ. Severe (weak) years are also associated with high (low) values of P_exc and, except for 26 May and 15 June 2007, the groups of severe and weak years are fully separated by contour P_exc=0.5 (which, in fact, corresponds to the probability of exceedance of μ_BA in the null model, where meteorological factors are not taken into account). Values of P_exc associated with severe years tend to gradually increase along the pre-fire season and, on 30 June, five out the six severe years present values above 0.7 (a threshold that is barely surpassed by just two of the 27 moderate years). Finally, it is worth noting that the contour lines of P_exc(ψ, χ) become steeper along the pre-fire season, a result in agreement with the steep increase in parameter a (Fig. 3).

Figure 4 Scatter plots of the 39-year sample in the space of covariates ψ(d) vs. χ for (a) 26 May (d=56), (b) 15 June (d=76) and (c) 30 June (d=91). Groups of severe (weak) years are identified by the red (green) circles. The straight lines are contours of P_exc(ψ, χ).

Results obtained suggest that the likelihood of a given year to belong to the severe or to the weak groups may be estimated based on the value of P_exc(d) of each year as estimated by the fitted model on chosen day d (Fig. 5). For a given year, values of P_exc(d) change very slowly from day to day. This is to be expected since (1) the models are all fitted to the same sample of log ₁₀BA, (2) covariate χ is the same in all models and (3) covariate ψ(d) has high serial correlation since, according to Eq. (1), index D_pfs(d) accumulates values of DSR since 1 April up to the considered day. Therefore, any decision taken for a given year, based on the respective estimate of P_exc(d) from the model fitted on day d, is not expected to drastically change in the next few days unless there is an incoming sequence of very high (or very low) daily values of DSR that will considerably change ψ(d) and therefore P_exc(d).

Figure 5 Temporal evolution of P_exc(ψ, χ) during the pre-fire season (26 May to 30 June). Severe (weak) years are identified by the red (green) curves. The black horizontal lines represent thresholds used in Type A (P_exc=0.5) and Type B (P_exc=0.7) decisions. Vertical dashed lines delimit the phase where Type A and Type B decisions were checked.

The following two types of decisions are tested.

• Type A. If P_exc>0.5 at day d then the year is not classified as weak (i.e. it is either moderate or severe); otherwise, if P_exc≤0.5 then the year is not classified as severe (i.e. it is either moderate or weak).

• Type B. If P_exc>0.7 at day d then the year is classified as severe.

Results obtained indicate that decision types A and B should be applied to different periods of the pre-fire season from 26 May to 30 June (Table 1) and from 18 to 30 June (Table 2). Furthermore, for each type of decision, two phases are identified in the respective periods, the days of each phase being characterized by the same decisions taken for each year. In the case of Type A decisions, there is just one wrong decision (2007 is incorrectly classified as not being a weak year) during phase A1 (26 May to 22 June) and all decisions are correct during phase A2 (23 to 30 June). In the case of Type B decisions, and for the entire period (18 to 30 June), all severe years but one (1998) are correctly classified as severe; conversely, three moderate years are incorrectly classified as severe during phase B1 (18 to 22 June) and this decreases to two during phase B2 (23 to 30 June).

Table 1 Performance assessment of Type A decisions during the pre-fire season (26 May to 30 June) based on the model with two covariates.

Table 2 As in Table 1 but for performance assessment of Type B decisions for 18 to 30 June.

3.3 Model with one covariate

Despite its usefulness in characterizing the role played by meteorological factors during the pre-fire and the fire seasons, the model discussed in the last subsection cannot be used to anticipate the likelihood of a given fire season given that one of the two covariates (χ) is derived from daily information during that same fire season. However, results obtained in the previous subsection indicate that the role played by covariate ψ(d) in each fitted model becomes more and more relevant when day d progresses along the pre-fire season as suggested by the steady increase in coefficient a(d) that even becomes larger than b(d) after 25 June (Fig. 3), as well as by the increase in slope of contour lines of P_exc(d) (Fig. 4).

Therefore, for each considered day d, we fitted to the sample of log ₁₀BA the following normal model, now with the mean linearly depending just on covariate ψ(d):

$$\begin{array}{}\text{(13)}& {\log }_{\mathrm{10}}\mathrm{BA}\sim N({\log }_{\mathrm{10}}\mathrm{BA};p\left(d\right)\times \mathit{\psi }\left(d\right)+q,{\mathit{\sigma }}_{\mathrm{BA}}).\end{array}$$

Maximum likelihood estimates of model parameters p(d) and q are obtained using Eqs. (9) and (10) and it may be noted that again parameter q is the same for all fitted models, with q=μ_BA because ψ(d) has a zero mean. As in the case of coefficient a in the model with two covariates, coefficient p(d) increases as the considered day d progresses along the pre-fire season (Fig. 6). For each day d the relative role played by covariates ψ(d) and χ may be assessed by comparing the log-likelihood ratio statistics $-\mathrm{2}\ln \left(L_{\mathrm{0}}/L_{\mathrm{1}}\right)$ and $-\mathrm{2}\ln \left(L_{\mathrm{1}}/L_{\mathrm{2}}\right)$, where L₀, L₁ and L₂ are the likelihood functions of the null model, the model with one covariate and the model with two covariates (Fig. 7). These two ratios represent the increases in likelihood of the sample of BA when replacing the null model (with no covariates) by the model with covariate ψ and then the latter model by the model with covariates ψ and χ.

Figure 6 Temporal evolution of coefficient p of the normal model with one covariate.

Values of $-\mathrm{2}\ln \left(L_{\mathrm{0}}/L_{\mathrm{1}}\right)$ increase as the considered day d progresses along the pre-fire season, contrasting with the behaviour of $-\mathrm{2}\ln \left(L_{\mathrm{1}}/L_{\mathrm{2}}\right)$ where a decrease (albeit more moderate) is observed. The former ratio even becomes larger than the latter from 23 to 30 June; however, it may be noted that, for models fitted on days between 26 May and 13 June, values of $-\mathrm{2}\ln \left(L_{\mathrm{0}}/L_{\mathrm{1}}\right)$ are smaller than the critical value at the 5 % significance level, indicating that the null model is to be retained (against the alternate model with one covariate). Nevertheless, the model with one covariate is still tested along the entire pre-fire season (26 May to 30 June) because, as shown in Fig. 4a, b, the severe (weak) groups tend to present for the most part high (low) values of covariate ψ. Reinforcing the relevance of this covariate when estimating the likelihood of a given year to belong to the severe class, we find that, when, for each considered fixed day d, we restrict to years with ψ(d) larger than the respective daily median, there is a positive correlation between ψ(d) and χ that, except from 5 to 12 June, is significant at the 5 % level (Fig. 8).

Figure 7 Temporal evolution of the log-likelihood ratio statistics $-\mathrm{2}\ln \left(L_{\mathrm{0}}/L_{\mathrm{1}}\right)$ between the null model and the model with one covariate and $-\mathrm{2}\ln \left(L_{\mathrm{1}}/L_{\mathrm{2}}\right)$ between the model with one covariate and the model with two covariates. The horizontal dashed line represents the critical value for the statistic at the 5 % significance level.

Figure 8 Temporal evolution of the correlation between ψ(d) and χ from 26 May to 30 June for all years (dotted curve) and restricting to years ψ(d) larger than the respective median (solid curve). Black circles in the solid curve identify values of correlation that are significant at the 5 % level.

As with the model with two covariates, a given year is anticipated as belonging to the severe or the weak groups by making Type A and Type B decisions based on the value of P_exc(d) of each year as estimated by the fitted model (with one covariate) on that fixed day d. It is worth noting the use of the wording “is anticipated as” (instead of “is classified as” employed in the previous section), which is meant to enhance the prognostic character of the model with covariate ψ. Since the model with one covariate has lower variability than the model with two covariates, the threshold of 0.7 (previously used in Type B decisions) is now lowered to 0.66 (Fig. 9).

Figure 9 Temporal evolution of P_exc(ψ) during the pre-fire season (26 May to 30 June). Severe (weak) years are identified by the red (green) curves. The black horizontal lines represent thresholds used in Type A (P_exc=0.5) and Type B (P_exc=0.66) decisions. Vertical dashed lines delimit the phase where Type A and Type B decisions were checked.

Three phases (A1, A2 and A3) are used to assess Type A decisions (Table 3). In all phases, only 1998 (out of the six severe years) is incorrectly anticipated as a non-severe year. In turn, the number of weak years incorrectly anticipated as non-weak decreases from three in phase A1 (26 to 30 May) to two in phase A2 (31 May to 5 June) and then to just one in phase A3 (6 to 30 June). Assessment of Type B decisions (Table 4) is performed in two phases (B1 and B2). Regarding the severe years correctly classified as severe, it is worth pointing out that, apart from phase B1 beginning 3 d later (on 21 June instead of 18 June), there is no decrease in performance from the model with two covariates, 1998 being again the only missed severe year in phases B1 and B2. Finally, no virtual decrease is found in performance in the number of moderate years incorrectly anticipated as severe, that, as with the model with two covariates, decrease from three in phase B1 (21 to 26 June) to two in phase B2 (27 to 30 June).

Table 3 As in Table 1 but for the model with one covariate.

Table 4 As in Table 2 but for the model with one covariate.

4 Discussion and conclusions

The increasing dimension of the impacts of extremely large wildfires that have been affecting Portugal in the last decades calls for the improvement of diagnostic and prognostic tools designed to assist forest and fire managers in making better decisions on both prevention and combat. For this purpose, we set up a model that allows outlooks, up to 1 month ahead, of wildfire potential in Portugal during the fire season (defined as July and August).

The model consists of a normal distribution of the decimal logarithm of the yearly cumulated values of BA during the fire season (log ₁₀BA) where the mean of the distribution linearly depends on a covariate ψ(d), defined at a given fixed day d of the pre-fire season, which responds to the accumulation of thermal and water stress up to the chosen day d. This model with one covariate results from the simplification of a model with two covariates ψ(d) and χ, the latter covariate being sensitive to the occurrence of hot and dry spells associated with strong winds that are the key triggering mechanism for the onset and spread of very large fire events (Amraoui et al., 2013). The rationale behind setting up the model with covariate ψ(d) is twofold: (1) long periods of warm and dry weather during the pre-fire season induce high levels of vegetation stress that increase the probability of occurrence of large fire events as shown for Portugal (Trigo et al., 2006) and even the Mediterranean basin (Gudmundsson et al., 2014; Turco et al., 2017); (2) in fact, it is now well-known that soil moisture deficit and drought (associated with large values of ψ) have an impact on hot summer extremes in Mediterranean regions (Vautard et al., 2007; Hirschi et al., 2011), a feature that fostered the development of both statistical and dynamical seasonal fire forecasting approaches in the Mediterranean region (Gudmundsson et al., 2014; Turco et al., 2017, 2018).

Performance of the model was evaluated for the 39-year period 1980–2018 by testing two types of decisions based on the temporal evolution of P_exc: for Type A, the decision is on whether a given year is classified as not being severe or as not being weak, whereas for Type B it is on whether the year is classified as severe. Type A decisions are made during the entire pre-fire season (26 May to 30 June) and Type B decisions are restricted to the end of the pre-fire season (21 to 30 June) when the role played by the cumulation of vegetation stress becomes as relevant as the one by the occurrence of extreme weather events during the fire season. In the case of Type A decisions, only 1998 (out of the six severe years) was incorrectly anticipated as non-severe. For Type B decisions, all severe years but 1998 were correctly anticipated as severe, reflecting the prominent role of cumulated vegetation stress at the end of the pre-fire season in favouring the occurrence of severe years. However, a test similar to Type B, but applied to weak years, has poor performance and this may be viewed as an indication that low values of ψ along the pre-fire season have a moderate impact on the likelihood of a following fire season with a low value of χ and therefore on the likelihood of having a weak year.

Since both Type A and Type B decisions have shown to be incorrect for the fire season of 1998, the temporal evolution of daily values of DSR is worth being analysed along the pre-fire and the fire seasons (Fig. 10). Except for two peaks in the last 10 d of June, the pre-fire season of 1998 is dominated by moderate daily values of DSR which are very close to zero between 10 May and early June, a period characterized by significant amounts of rainfall. Therefore, by the end of May, vegetation was not subject to either thermal or water stress and this translates into low values of ψ(d) along the pre-fire season (Fig. 4). This unstressed condition of vegetation explains why the two large peaks of daily DSR in the second half of June did not trigger any large fire events as well as why the large peak of DSR at the end of the first half of July was followed by days of moderate burned area. However, the meteorological conditions drastically changed in August, with the first half of the month being dominated by values of daily DSR well above the median that culminated with a sequence of days with DSR much higher than percentile 90, and the second half presenting a peak of DSR very close to percentile 90. The impact of this high number of extreme days translates into a large value of χ for 1998, the seventh largest in the 39 years analysed and the third of the six severe years. Striking a vegetation stressed by the warm and dry conditions that prevailed since the beginning of August, the two peaks of DSR triggered two sequences of days affected by large fire events that made 1998 rank fifth in cumulated BA in 1980–2018.

Figure 10 Temporal evolution of daily values of DSR for 1998 (red curve) and of the daily values of the median (solid black curve) and of percentile 90 (dotted black curve) for the period 1980–2018. The dark grey bars indicate the cumulated daily values of BA (in hectares) for Portugal. The vertical black line subdivides the period into pre-fire and fire seasons.

As shown in Fig. 4, 1998 is the only severe year when a value of ψ at the last day of the pre-fire season (30 June) below the mean minus 1 standard deviation is followed by a value of χ above the mean plus 1 standard deviation. The exceptional character of 1998 is reinforced by the fact that, even when considering all 39 years, this condition is fulfilled in only 2 years, i.e. 1998 and 1985 (a moderate year). However, despite the extreme meteorological conditions observed in August, 1998 just ranks fifth among the six severe years, suggesting that the mild conditions during the pre-fire season still had a mitigating role. Conversely, the exceptionality of 1998 puts an emphasis on the importance of short-term extreme events as triggers of large-fire events, a feature that is not considered when restricting to the information provided by maps of long-term seasonal forecasts of temperature and rainfall anomalies, undermining the usefulness of current seasonal forecasting outputs for wildland fire potential outlooks.

Several improvements to the proposed prognostic model (with covariate ψ) are currently being considered. First, no advantage was taken of the fact that covariates ψ and χ are positively correlated (when ψ is larger than the median). This information may be added by setting up a model of the distribution of χ using information about meteorological conditions along the pre-fire season and then incorporating this information in the prognostic model. Second, the model was developed for the entire territory of Portugal, not taking into account regional characteristics of climate, land cover and fire regime. It is therefore worth setting up regional models over areas with distinct pyrogeographical characteristics. Finally, the outlooks with respect to the months of July and August and given the observed trend to have a longer fire season it would be very useful to extend them until October. For instance, as pointed out by Beighley and Hyde (2018), 36 % of the total area burned in 2009–2017 was outside the period from July to September, a fraction that is 3 times larger than the value of 12 % observed when considering the period 2001–2008.

A beta version of the proposed model has been experimentally running since May 2017 and the correctness of the outlooks (made in real time) about 2017 being a severe year (based on a Type B decision) and about 2018 not being severe (based on a Type A decision) give confidence about the potential of the model to be operationally used to anticipate the occurrence of severe years.

Daily outlooks are currently available at CeaseFire (http://idlcc.fc.ul.pt/CeaseFire/index.php, last access: 15 July 2019), a website designed to integrate and disseminate relevant meteorological information to the user fire community by means of a simple, fast and user-friendly interface (Evans, 2018). Developed by Instituto Dom Luiz (IDL) at the Faculty of Sciences of the University of Lisbon (Portugal), the platform relies on data provided in near-real time by LSA SAF, the EUMETSAT Satellite Application Facility for Land Surface Analysis (Trigo et al., 2011). Currently there are about 900 registered users in the CeaseFire platform, most of them from national authorities and services, as well as from municipalities and private companies, namely from the paper and pulp industry.