The Signature of Climate in Annual Burned Area in Portugal

Abstract

Portugal is by far the country most affected by wildfires in Mediterranean Europe. The increase in frequency and severity of extreme years in the last two decades calls for a better understanding of the role played by climate variability and climate change. Using data covering a period of 44 years (1980–2023), it is shown that the distribution of annual burned area in Portugal follows a Rayleigh distribution whose logarithm of the scale parameter depends linearly on Cumulative Daily Severity Rate ($CDSR$). Defined for each year as the sum of the mean Daily Severity Rate over Portugal from 1 January to 31 December, $CDSR$ is a measure of the dryness of dead fuels as induced by atmospheric conditions. Changes along the years of the modeled average explain 56% of the interannual variability of the annual burned area. When comparing the model for 30-year subperiods 1980–2009 and 1994–2023, large decreases are observed in return periods of annual burned area amounts, from 35% for amounts greater than 120 thousand hectares up to 49% for amounts greater than 200 thousand hectares. The proposed model is a useful tool for fire management under present and future climate conditions.

1. Introduction

From 1980 to 2023, wildfires in mainland Portugal burned a total of 4,987,678 hectares [1], an area equivalent to 54% of the territory. During that 44-year period, Portugal was recurrently affected by extreme fire seasons of increased severity. Thresholds of 180, 470, and 530 thousand hectares of the annual burned area were successively broken in the decades of 1990, 2000, and 2010, respectively, in 1991 (182,486 ha), 2003 (471,750 ha), and 2017 (539,921 ha), with 2003 being closely followed, in 2005, by another extreme year (346,718 ha). The amount of total burned area and the severity of the extreme fire seasons make Portugal by far the country most affected by wildfires among the five southern countries of Europe (France, Greece, Italy, Portugal, and Spain).

Wildfire activity in Portugal is mainly associated with large-scale atmospheric patterns [2,3,4,5] and synoptic-scale circulation weather types [6,7] that, by steering easterly winds carrying hot and dry air over the country, set up an appropriate background for the spread and buildup of intense fires [8,9,10]. Wildfire activity is also associated with vegetation stress as induced by prolonged warm and dry conditions that decrease the moisture content of deep, compact organic layers [11,12], as well as by heat waves that dry the moderate duff layers and medium-size woody material and turn flammable the litter and the fine fuels [13]. The role played by drought and heat spells is especially prominent in extreme years; the drought of 2005 was one of the worst episodes ever recorded [14], and the fire seasons of 1991 and 2003 were punctuated by very strong heat waves; the spell of 8–22 July 1991 has affected the whole territory [15] and that of 1–15 August 2003 consisting of an exceptionally intense event [16]. During the record year of 2017, Portugal was affected by a very severe drought and by two strong heatwaves, in 7–24 June and 12–17 July, respectively. As in other regions of the world [17,18], the compound drought and heat wave events strongly contributed to the exceptionality of 2017, both in terms of amount, which represented almost 60% of the total burned area in Europe in that year, and in terms of the extent of the fire season, which was marked by two tragic events, the first one in 17–20 June and the second one in 15–17 October [19,20]. With 53 thousand hectares of burned area and a toll of 66 dead and 254 wounded, the event of mid-June was the deadliest ever recorded in Portugal; in turn, during the event of mid-October that accounted for 50 dead and 71 wounded, the burned area reached the all-time record of 220 thousand hectares in 24 h, an amount that is twice the annual average for the period 1980–2023. Both events were exacerbated by a combination of extremes; in the case of the mid-June event, the effects of anomalous high values of air temperature and low values of air relative humidity associated with the compound drought and heat spell were exacerbated by very unstable atmospheric conditions that steered downburst events [19]. The case of mid-October was the result of a combination of extreme meteorological conditions with negligent human behavior; the hydric vegetation stress induced by the prolonged drought conditions was magnified by the intense southerly winds associated with the passage of hurricane Ophelia off-coast of Portugal, and there was a record number of ignitions set up by farmers in a critical day of fire danger [20].

Wildfire activity in Portugal is to be viewed as the result of a complex interplay among climate, landscape, and human activity factors, and therefore accounting for the whole process requires sound information about weather patterns [21,22], vegetation dynamics [11,23], fire activity [24], and burned area [25]. Remote sensing is playing an increasing role in wildfire activity monitoring, namely the retrieval of key parameters like land surface temperature [26,27,28,29,30] and emissivity [31,32] and the estimation of surface radiative balance [33], as well as hot spot detection [15,34], burned area mapping and dating [35,36,37,38,39], and post-fire vegetation recovery [40,41,42].

The increasingly adverse impacts of wildfires at the economic, social, environmental, and ecological levels have boosted the development of products tailored to the needs of fire managers, such as the set of fire products [15,43] developed by the Land Surface Analysis Satellite Application Facility (LSA SAF) that is part of the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Application Ground Segment [44]. Special attention has been devoted to the development of 1- to 10-day forecasts of fire danger [45,46,47] and to seasonal outlooks of wildfire potential [48,49,50] that are disseminated on a daily basis by the LSA SAF and the University of Lisbon [51,52]. Fire danger is most often assessed based on indices derived from meteorological parameters like those that compose the Canadian Forest Fire Weather Index System (CFFWIS), which is especially appropriate to the Mediterranean ecosystems [53] and has been adopted by the Portuguese Weather Service (IPMA) since 1998 [54]. The Fire Weather Index ($FWI$), the last component of CFFWIS, rates fire intensity and has been extensively used to define classes of fire danger [45]. Directly derived from $FWI$, the Daily Severity Rating ($DSR$) grades the difficulty of controlling fires and has proven to be useful for outlooks of the severity of the fire season [48].

The issue of climate change, whose impacts are especially adverse for Mediterranean ecosystems [55], has called for studies dedicated to the role played by climate on the current [56,57] and future [58,59] dynamics of the fire regimes. The aim of the present study is to assess the impacts of climate variability and climate change on the annual burned area in Portugal, considering a period of almost four decades and a half (1980–2023). For this purpose, it is first shown that the burned area follows a Rayleigh distribution with Cumulative Daily Severity Rate ($CDSR$) as a covariate of the scale parameter of the distribution. $CDSR$ is defined on a yearly basis as the accumulated value (from 1 January to 31 December) of the daily spatial averages of $DSR$ over Portugal. The impact of climate on the annual burned area is then evaluated based on the temporal evolution of the distribution of burned area associated with the variability of the scale parameter.

A brief description of the data used in this study is provided in the next section, together with the key aspects of the theoretical background related to the Rayleigh distribution and its fitting to the data, tests of goodness of fit, and trend analysis. The results obtained are presented in Section 3 and then discussed in Section 4. Final considerations, at a more general level, are provided in Section 5.

This paper is the outcome of the “2023 Academic Career Award in Meteorology” that the author received from the Portuguese Association of Meteorology and Geophysics (APMG). In line with this award, the topic of the paper was chosen so that a review could also be made of the author’s rather long academic career that he started on the 1 August 1982 as a teaching assistant at the Department of Physics of the Faculty of Sciences of the University of Lisbon, Portugal. This is the only reason for the unusually long list of co-authored papers that are referred to in the Introduction and Concluding remarks sections. The list is just to be viewed as a way of paying tribute to all those who have contributed to making up the author’s career.

2. Data and Methods

Data used in this study cover the period 1980–2023 and consist of annual amounts of burned area in Portugal (hereafter $BA$) and daily values of $DSR$.

The time series of $BA$ is from the official database of the Portuguese Institute for Nature Conservation and Forests (ICNF) and is available in the reports of the European Forest Fire Information System (EFFIS) [1]. Statistical information about wildfire activity at the regional and national levels together with records of individual fire events for the period 1980–2023 are also available at the ICNF site [60]. Under Portuguese law, firefighters are the source of information about each rural fire occurring in Portugal. During the period 1980–2000, location and burned area verification was done by the Portuguese “Forest Guards”, and, after 2001, the verification procedure was improved by introducing SGIF, the Information Management System for Forest Fires [61].

Daily values of $DSR$ were extracted from the Copernicus Climate Change Service and are available online [62]. $DSR$ is directly obtained from $FWI$ [48] according to the relationship:

$$DSR=0.0272FW{I}^{1.77}.$$

(1)

By weighting $FWI$ sharply as it increases, transformation (1) makes $DSR$ more suitable than $FWI$ to be averaged in space and accumulated/averaged in time; for instance, monthly averages and accumulated values since 1 April of $DSR$ averaged over Portugal have been successfully used to develop statistical models of the logarithm of burned area amounts accumulated in the fire season [48,58]. Daily values of DSR are defined in a $0.5°\times 0.5°$ regular lat-lon grid, with 150 grid points covering the Portuguese territory. In this study, for each day of the period 1980–2023, averages were computed of $DSR$ on grid points over Portugal. Then, for each year, the Cumulated Daily Severity Rate (CDSR) was obtained by adding the daily spatial averages of DSR from 1 January to 31 December. For each year, daily spatial averages of $DSR$ over Portugal were also accumulated from 1 January up to each day of the year; the notation $CDSR\left(d\right)$ will be used to refer to the cumulative $DSR$ from 1 January to day $d$.

The statistical models in this study are based on the Rayleigh distribution [63], whose support is $[0,\mathrm{\infty })$ and is, therefore, appropriate to model non-negative variables such as $BA$. The probability distribution function (pdf) is as follows:

$$f\left(BA;\sigma \right)={\displaystyle \frac{BA}{{\sigma }^2}}{e}^{{-\left(BA\right)}^2/\left(2{\sigma }^2\right)},$$

(2)

where $\sigma$ is the scale parameter of the distribution.

The cumulative distribution function (cdf) is as follows:

$$F\left(BA;\sigma \right)=1-e^{-{\left(BA\right)}^2/\left(2{\sigma }^2\right)}.$$

(3)

The mean, ${BA}_{mean}$, the median, ${BA}_{median}$, and the mode, ${BA}_{mode}$, of the Rayleigh distribution, respectively, are as follows:

$${BA}_{mean}=\sigma \sqrt{\pi /2}\simeq 1.25 \sigma$$

(4a)

$${BA}_{median}=\sigma \sqrt{2\ln 2}\simeq 1.18 \sigma$$

(4b)

$${BA}_{mode}=\sigma .$$

(4c)

Since ${BA}_{mode}<{BA}_{median}<{BA}_{mean}$, the distribution is skewed to the right, and the mean to median ratio is as follows:

$$MMR={\displaystyle \frac{{BA}_{mean}}{B{A}_{median}}}={\displaystyle \frac{1}{2}}\sqrt{\pi /\ln 2}\simeq 1.06.$$

(5)

The standard deviation, ${BA}_{std}$, of the Rayleigh distribution is as follows:

$${BA}_{std}=\sigma \sqrt{(4-\pi )/2}\simeq 0.66 \sigma$$

(6)

and percentile q, ${BA}_q$, is as follows:

$${BA}_q=\sigma \sqrt{-2\ln (1-q/100)}.$$

(7)

As expected, Equation (4b) follows from Equation (7) by setting $q=50$; on the other hand, it follows from Equations (4a) and (7) that:

$${BA}_q={BA}_{mean}\Rightarrow -2\ln \left(1-q/100\right)=\pi /2\Rightarrow q/100=1-e^{-\pi /4}\simeq 0.54,$$

(8)

i.e., the mean of the Rayleigh distribution virtually coincides with percentile $54$ (${\mathrm{i}.\mathrm{e}., B{A}_{mean}\cong BA}_{54}$).

Considering Equations (4a) and (6), the coefficient of variation is as follows:

$$CV={\displaystyle \frac{{BA}_{std}}{B{A}_{mean}}}=\sqrt{(4-\pi )/\pi }\simeq 0.52.$$

(9)

The return period of the event ${BA=BA}^{*}$ is as follows:

$$RP\left({BA}^{*}\right)={\displaystyle \frac{1}{1-F\left({BA}^{*};\sigma \right)}}=e^{{-\left({BA}^{*}\right)}^2/\left(2{\sigma }^2\right)}.$$

(10)

Given a sample $\left\{{BA}_1\cdots {BA}_N\right\}$ of $N$ events, the maximum likelihood estimate of $\sigma$ is given by the following:

$$\widehat{\sigma }=\sqrt{{\displaystyle \frac{1}{2N}}{\sum }_{i=1}^N{\left(B{A}_i\right)}^2},$$

(11)

Quality of fit is confirmed by examining q-q plots of empirical versus modeled quantiles, and goodness of fit is checked using the Anderson-Darling test to assess the strength of the evidence against the null hypothesis that the sample follows a Rayleigh distribution; confidence levels are obtained by randomly generating 1000 data samples from the Rayleigh distribution with estimated $\widehat{\sigma }$ [45].

In a previous study for Portugal [48], it is shown that the distribution of the logarithm of the annual burned area in July and August follows a normal distribution whose mean linearly depends on accumulated values of $DSR$ since 1 April. This suggests incorporating $CDSR$ as a covariate of the scale parameter of the Rayleigh distribution according to the linear relationship:

$$\ln \sigma =m\times CDSR+b.$$

(12)

Applying logarithms to both members of Equations (4a) and (7) and taking Equation (12) into account, the following expressions are obtained:

$$\ln {BA}_{mean}=m\times CDSR+\left(b+{\displaystyle \frac{1}{2}}\ln (\pi /2)\right)$$

(13a)

$$\ln {BA}_q=m\times CDSR+\left(b+{\displaystyle \frac{1}{2}}\ln [-2\ln \left(1-q/100\right)]\right)$$

(13b)

meaning that when $\ln BA$ is plotted against $CDSR$, the dependence of the mean and any quantile on $CDSR$ is represented by parallel lines with slope $m$.

Let $\left\{C{DSR}_1\cdots {CDSR}_N\right\}$ be the values of $CDSR$ associated to the sample $\left\{{BA}_1\cdots {BA}_N\right\}$. Estimates $\widehat{m}$ and $\widehat{b}$, respectively of slope $m$ and intercept $b$ in Equation (12), may be obtained either from Equation (13a) by least squares regression or from Equation (13b) by quantile regression; confidence intervals for estimated parameters may be obtained from bootstrap surrogates.

In case $\left\{C{DSR}_1\cdots {CDSR}_N\right\}$ are independent, the associated pdf is as follows:

$$\overline{f}\left(BA;{\sigma }_1\cdots {\sigma }_N\right)={\displaystyle \frac{1}{N}}{\sum }_{i=1}^{N}f(BA;{\sigma }_i)={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{\displaystyle \frac{BA}{{\sigma }_i^2}}{e}^{{-\left(BA\right)}^2/\left(2{\sigma }_i^2\right)},$$

(14)

where ${\sigma }_i=\mathrm{e}\mathrm{x}\mathrm{p}(m\times CDS{R}_i+b)$; and the cdf is as follows:

$$\overline{F}\left(BA;{\sigma }_1\cdots {\sigma }_N\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^{N}F(BA;{\sigma }_i)={\displaystyle \frac{1}{N}}\sum _{i=1}^N\left[1-e^{-{\left(BA\right)}^2/\left(2{\sigma }_i^2\right)}\right].$$

(15)

Goodness of fit of distribution $\overline{f}\left(BA;{\sigma }_1\cdots {\sigma }_N\right)$ is again assessed based on q-q plots and using the Anderson–Darling test.

The existence of trends in time series is assessed using the Mann–Kendall test of the null hypothesis of trend absence against the alternative of existence of trend; in case the null hypothesis is rejected, the trend is estimated by performing Theil–Sen robust linear regression on data [64].

The statistical model of $BA$ is derived in two steps. First, it is shown that, without the three largest values, the sample of $BA$ covering the period 1980–2023 follows a Rayleigh distribution with constant scale parameter $\sigma$; this is done by showing that the estimates of $MMR$ and $CV$ of the sample are very close to the theoretical values for a Rayleigh distribution as given by Equations (5) and (9). In the second step, it is shown that the entire sample follows a Rayleigh distribution whose logarithm of $\sigma$ linearly depends on $CDSR$. This is done by using quantile regression to show that, in accordance with Equation (13a), percentiles 25 and 75 of the logarithm of $BA$ depend linearly on $CDSR$ with virtually the same slope; then, using quantile and least squares regressions, it is shown that, in accordance with Equations (8) and (13b), percentile 54 and the mean of the logarithm of $BA$ linearly depend on $CDSR$, with the average of the respective slopes virtually coinciding with the common slope estimated for percentiles 25 and 75.

The fitted model is finally used to assess the impacts of climate variability and climate change on $BA$. The impact of climate variability is characterized by computing the contribution of the time series of the annual means of the fitted model to the total variance of the time series of $BA$. The impact of climate change is in turn assessed by comparing the return period of several thresholds of annual burned area for two subperiods of 1980–2023, covering respectively the first thirty years (1980–2009) and the last thirty years (1994–2023).

3. Results

3.1. BA and CDSR

Time series of $BA$ in 1980–2023 (Figure 1a) presents high interannual variability, ranging from 19 897 ha in 2008 up to 539 921 ha in 2017. As shown in

Table 1. Measures of location and dispersion respecting $BA$ in 1980–2023 and in subperiods 1980–2001 and 2002–2023.

Title 1	1980–2023 *	1980–2001	2002–2023
Mean ($B{A}_{mean}$)	113,356 [88,519] ha	92,450 ha	134,263
Median ($B{A}_{median}$)	86,287 [77,323] ha	83,095 ha	87,916
1st quartile ($B{A}_{25}$)	46,195 [44,496] ha	49,963 ha	42,084 ha
2nd quartile ($B{A}_{75}$)	143,604 [132,450] ha	137,252 ha	151,370 ha
Standard deviation ($B{A}_{std}$)	106,343 [49,023] ha	49,314 ha	140,743
Mean to median ratio ($MMR$)	1.31 [1.14]	1.11	1.53
Coefficient of variation ($CV$)	0.94 [0.55]	0.53	1.05

* Values in brackets are for 1980–2023 without the extreme years (2017, 2003, and 2005).

, the sample of $BA$ is skewed to the right, the mean being considerably larger than the median as indicated by the mean-to-median ratio ($MMR=1.31$). Reflecting the high interannual variability of the time series of $BA$, the coefficient of variation is close to unity ($CV=0.94$), meaning that the standard deviation is close to the mean.

The three largest values of $BA$ (in 2017, 2003, and 2005) are outstandingly larger than the remaining values, and this feature is especially conspicuous when representing the amounts of $BA$ in ascending order (Figure 1b). The considerable impact of the three extreme years on the statistics of $BA$ is evident when comparing the measures of location and dispersion of $BA$ in 1980–2023, with and without the extreme years (Table 1). According to the Mann-Kendall test, the null hypothesis of the absence of a trend in the time series of $BA$ cannot be rejected at the 5% significance level (p-value = 0.9).

Despite the absence of a trend, the first (1980–2001) and the second (2002–2023) halves of the period present quite different behavior, namely in terms of variability (Figure 1a), and this translates into changes in the measures of both location and dispersion (Table 1). There is a substantial increase in the mean, median, and percentile 75 from the first to the second half of the period, which indicates a shifting of the $BA$ distribution towards larger values. In turn, there is a sharp decrease in percentile 25 that, together with the increase of percentile 75, indicates an increase in dispersion, and this result is in line with the sharp increases in standard deviation and $CV$. There is also a considerable increase in skewness that translates into a substantial increase in $MMR$.

The time series of $CDSR$ in 1980–2023 (Figure 2a) presents a positive trend with large fluctuations. The positive trend is statistically significant at the 5% level (p-value = 0.002) according to the Mann–Kendall test and the trend slope is $14.4/\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$, as estimated using the non-parametric Theil–Sen robust linear regression. Values of $CDSR$ range from 1407 in 1997 to 3213 in 2017; the dry years of 1981, 1991, 1995, 2005, 2012, 2017, and 2023 being associated to peaks of $CDSR$. The two top values of $CDSR$, which coincide with the extreme years of 2017 and 2005, are outstandingly large, this feature becoming clear when the values of $CDSR$ are represented in ascending order (Figure 2b). The time series of $CDSR$, like that of $BA$, presents a larger variability in the second half of 1980–2023 than in the first half (Figure 2a).

3.2. Statistical Model of BA

When the extreme years of 2017, 2003, and 2005 are excluded from the $BA$ dataset, the obtained estimates of 1.14 for $MMR$ and 0.55 for $CV$ (Table 1) are close to the corresponding values of 1.06 and 0.52 for the Rayleigh distribution, shown in Equations (5) and (9).

It is, therefore, reasonable to consider the Rayleigh distribution as an appropriate model of $BA$ when excluding the extreme years. The obtained maximum likelihood estimate for the scale parameter is $\widehat{\sigma }=\mathrm{71,345} \mathrm{ha}$, and plots of the pdf and cdf of the fitted Rayleigh distribution are shown in Figure 3a and Figure 3b, respectively. The q-q plot of the empirical quantiles of the $BA$ dataset versus the estimated quantiles by the Rayleigh model (Figure 3c) indicates a good fit of the model to the data, with most of the points lying close to the 1:1 line. The goodness of fit is further confirmed by the Anderson–Darling test performed on the $BA$ dataset, whose result indicates that the null hypothesis that the sample follows the Rayleigh distribution cannot be rejected at the 5% level (p-value = 0.6).

However, when considering the entire $BA$ dataset (1980–2023), the obtained estimates of 1.31 for $MMR$ and 0.94 for $CV$ (Table 1) considerably depart from the corresponding values for the Rayleigh distribution suggesting that the model must be improved by allowing the scale parameter $\sigma$ to depend on a suitable covariate. The fact that the two topmost values of $CDSR$ correspond to two extreme years (2017 and 2005) and that the values of $CDSR$ for the three weakest years are among the eight lowest (Figure 2b) makes of $CDSR$ a good candidate for covariate of $\sigma$.

Following the procedure described at the end of Section 2, the dependence of $\sigma$ on $CDSR$ was derived from estimates of the linear relationship between $CDSR$ and the logarithms of selected quantiles of $BA$ as well as between $CDSR$ and the logarithm of the mean of $BA$.

First, estimates of the linear dependence on $CDSR$ of the logarithm of first and third quartiles of $BA$ (i.e., $B{A}_{25}$ and $B{A}_{75}$) were obtained by quantile regression:

$$\ln {\widehat{BA}}_{25}=0.0011\times CDSR+8.56$$

(16)

$$\ln {\widehat{BA}}_{75}=0.0012\times CDSR+9.42,$$

(17)

and it is worth noting that the two lines are almost parallel (Figure 4a).

Estimates of the linear dependence on $CDSR$ of the logarithm of percentile 54 of $BA$ and the logarithm of the mean of $BA$ (i.e., $B{A}_{54}$ and $B{A}_{mean}$) were then obtained by quantile regression and by least squares regression, respectively:

$$\ln {\widetilde{BA}}_{54}=0.0013\times CDSR+8.86$$

(18a)

$$\ln {\widetilde{BA}}_{mean}=0.0010\times CDSR+9.17,$$

(18b)

As shown in Equation (8), in the case of the Rayleigh distribution $B{A}_{54}\cong B{A}_{mean}$, and obtained differences in slope and intercept between Equations (18a) and (18b) reflect the different sensitivity of quantile and least square regressions; the common dependence of $B{A}_{54}$ and $B{A}_{mean}$ on $CDSR$ is therefore estimated by the bisector of lines defined by Equations (18a) and (18b):

$$\begin{gathered} \ln {\widehat{\mathit{BA}}}_{\mathit{mean}}=\ln {\widehat{BA}}_{54}=bisector\left(\ln {\widetilde{BA}}_{mean},\ln {\widetilde{BA}}_{54}\right) \\ =0.0012\times CDSR+9.02, \end{gathered}$$

(18c)

and again, it is worth noting that this new line and lines given by Equations (16) and (17) are almost parallel (Figure 4a).

The obtained linear dependence on $CDSR$ of the logarithm of the first and third quartiles and the mean of $BA$ as being described by almost parallel lines is consistent with the hypothesis of the distribution of $BA$ following a Rayleigh distribution whose logarithm of the scale parameter linearly depends on $CDSR$. Estimates of such dependence may be obtained from Equation (18c) using Equation (13a):

$${\ln \widehat{\sigma }}_{mean}=\ln {\widehat{\sigma }}_{54}=0.0012\times CDSR+8.79,$$

(19a)

as well as from Equation (16) using Equation (13b) with $q=25$:

$$\ln {\widehat{\sigma }}_{25}=0.0011\times CDSR+8.84,$$

(19b)

and from Equation (17) using Equation (13b) with $q=75$:

$$\ln {\widehat{\sigma }}_{75}=0.0012\times CDSR+8.91.$$

(19c)

The three estimates of the dependence of $\sigma$ on $CDSR$ are very close to each other (Figure 4b), and it is worth noting that all of them lie within the 95% confidence intervals for ${\widehat{\sigma }}_{25}$ and ${\widehat{\sigma }}_{75}$. Choosing the estimate given by Equation (19a), the dependence of the scale parameter of the Rayleigh distribution on $CDSR$ is given by the following:

$$\widehat{\sigma }=6572 {\mathrm{e}}^{0.0012\times CDSR}.$$

(20)

The pdf and cdf, $\overline{f}$ and $\overline{F}$, associated to the set of $CDSR$ values in 1980–2023 (Figure 5a,b) are obtained introducing Equation (20) into Equation (14) and Equation (15), respectively.

The goodness of fit is reflected in the q-q plot of the empirical quantiles of the $BA$ dataset versus the estimated quantiles of the Rayleigh model (Figure 5c), with all points except the extreme years (2017, 2003, and 2005) lying close to the 1:1 line. Moreover, the result of the Anderson–Darling test indicates that the null hypothesis that the sample follows the Rayleigh distribution with $CDSR$ as covariate cannot be rejected at the 5% level (p-value = $0.9$).

3.3. Interannual Variability of BA

The Rayleigh model with $CDSR$ as covariate allows computing the statistical distributions of $BA$ along the years (Figure 6) as well as evaluating the interannual variability of the different measures of location and dispersion. It is worth noting that, with the conspicuous exception of 2003, all values observed of $BA$ lie within percentiles 2.5 and 97.5 of the modeled distributions for the respective years, with 13 years corresponding to values below the first quartile, 18 years between the first and third quartiles, and 13 years (including 2003) above the third quantile. The time series of the mean explains 39% of the total variance of the time series of $BA$, and this amount increases to 57% when excluding 2003.

The ability of the model to capture the exceptional nature of the extreme years of 2017 and 2005 is worth emphasizing, with the recorded values of $BA$ corresponding to values of cumulative probability of 0.8 for 2017 and 0.7 for 2005 according to the respective statistical models. As already noted, this is not the case for 2003 since, according to the distribution modeled for that year, the cumulative probability for the recorded value is very close to 1.

Figure 7 presents, for the three extreme years (2017, 2005, and 2003), time-series of daily values of $CDSR\left(d\right)$ (Figure 7a) and $DSR$ (Figure 7b) for the period spanning from 1 June to 31 October; for each day, the range between daily percentiles 5 and 95 of $CDSR\left(d\right)$ and $DSR$ for the entire period of study (1980–2023) is also provided as a climatological background.

There is a sharp contrast between 2003 and the other two extreme years (2005 and 2017) for both $CDSR\left(d\right)$ and $DSR$. Contrary to 2005 and 2007, from 1 June to 31 October, present high values of $CDSR\left(d\right)$, above percentile 95 for virtually all days in 2005 and from October onwards in 2017, values of $CDSR\left(d\right)$ in 2003 are moderate throughout the entire period.

A contrast is also observed in the case of $DSR$ from June to September, with 2003 presenting several sequences of high peaks beyond percentile 95 (two located in June, two in August, and three in September) interleaved by long periods of low values of $DSR$ whereas 2005 and 2017 present peaks of lower magnitude and fluctuations of lower amplitude; in October, the three extreme years also show different behaviors, 2003 presenting low values of $DSR$ throughout the month, 2005 presenting high and low values in the first and second halves of the month, respectively, and 2017 presenting an extremely high peak of $DSR$ on 15 October that is preceded by a week with high values and followed by another week of high values at the end of the month.

The distinct character of 2003 compared with 2005 and 2017 is conspicuous when the respective values of $DSR$ are rearranged in descending order (Figure 8). Excepting for the outstandingly high value of $DSR$ on 15 October 2017, the first eight values of $DSR$ in 2003 overcome the corresponding ones of 2005 and 2017. On the other hand, when compared with 2005 and 2017, values of $DSR$ in 2003 increasingly deviate towards lower values, whereas those of 2005 and 2017 become closer and closer until rank 113. The different behavior of 2003 versus 2005 and 2017 is also reflected both in the mean values for the period (12.0 in 2003 versus 15.6 and 17.5 in 2005 and 2017) and in the standard deviation (8.8 in 2003 versus 8.2 and 7.4 in 2005 and 2017).

3.4. Trend in $\sigma$ and Return Periods of BA

Reflecting the temporal behavior of $CDSR$, the time series of scale parameter $\sigma$, as derived from Equation (20), present large fluctuations superimposed on a positive trend of $1148 \mathrm{h}\mathrm{a}/\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$ as estimated by a Theil-Sen robust linear regression (Figure 9a). The trend in $\sigma$ is statistically significant at the 5% level (p-value = 0.002) according to the Mann-Kendall test, and the two extreme years of 2005 and 2017 are associated with very large peaks of $\sigma$. The increase in both magnitude and variability of $\sigma$ along 1980–2023 translates into changes in the statistical distribution of $BA$ that are conspicuous when comparing the pdf $\overline{f}$ (Figure 9b) and the cdf $\overline{F}$ (Figure 9c) curves of the Rayleigh models associated to the sets of $CDSR$ in the first (1980–2009) and the last (1994–2023) 30 years of the period considered (1980–2023).

When going from the first to the last 30-year subperiods, although there is a small shift towards the right in the maxima of the two pdf curves (Figure 9b), corresponding to a moderate change in the respective modes from 52,100 ha to 58,500 ha (12% increase), the major changes are observed in the tails of the pdf curves that are considerably heavier in the last 30-year period, resulting in substantial increases in the means, from 88,679 ha to 117,080 ha (32% increase) and in the standard deviations, from 67,706 ha to 97,863 ha (45% increase). In the case of the cdf curves (Figure 9c), there is an overall shift to the right when going from the first to the last 30-year periods that leads to large increases in the median, from 72,900 ha to 90,800 ha (25% increase), as well as in the first and third quartiles, respectively, from 44,900 ha to 54,200 ha (21% increase) and from 112,200 ha to 147,000 ha (31% increase).

The displacement to the right of the cdf $\overline{F}$ curve for the last 30-year period compared to that for the first 30-year period indicates that, for a chosen threshold $B{A}^{*}$, the respective cdf values are lower in the last period than in the first; attending to relation (10) that defines return period, such decrease implies a decrease in the respective return periods for the chosen threshold. The decreases in return periods become larger in relative terms with increasing thresholds (Figure 9d), e.g., from 4.1 to 2.7 years (35% decrease) for amounts greater than 120,000 ha, from 8.1 to 4.3 years (47% decrease) for amounts greater than 160,000 ha, and from 19.8 to 8.1 years (59% decrease) for amounts greater than 200,000 ha.

4. Discussion

The annual amount of burned area in Portugal ($BA$) for the period 1980–2023 is characterized by very large interannual variability (Figure 1a). This feature is in good part due to the impact of three extreme years (2017, 2003, and 2005), whose exceptionality is clear when placing the yearly amounts in ascending order (Figure 1b).

Annual amounts of $CDSR$ for the same period of 1980–2023 also present large fluctuations about a positive trend (Figure 2a). Years affected by severe droughts are associated with peaks in the $CDSR$ time series, an indication that this quantity is appropriate to rate the long-term effects of meteorological conditions on the moisture content of organic layers that modulate the landscape proneness to be affected by large vegetation fires. For instance, the linkage between $BA$ and $CDSR$ translates into the fact that the two topmost values of $CDSR$ occur in the extreme years of 2005 and 2017 (Figure 2b), and years with large amounts of annual burned area (e.g., 1981, 1985, 1991, 1995, 2000, 2005, 2013, 2017, and 2022) are associated with large values of $CDSR$.

Although a Rayleigh distribution with constant $\sigma$ is adequate to model the distribution of $BA$ without the extreme years of 2017, 2003 and 2005, extending the model to the entire period of 1980–2023 requires incorporating $CDSR$ as a covariate of the scale parameter $\sigma$. The impact of $CDSR$ on $BA$ is accordingly assessed by modeling the distribution of $BA$ by means of a Rayleigh distribution with $CDSR$ as covariate of $\sigma$. Except for 2003, all annual amounts of $BA$ are within percentiles 2.5 and 97.5 of the distribution for the corresponding year.

The exceptionality of 2003 is conspicuous when analyzing daily values of both $CDSR\left(d\right)$ and $DSR$ (Figure 7 and Figure 8). When compared with the remaining years, the time series of $CDSR\left(d\right)$ along the extended summer (June–October) is characterized by moderate values well within the band delimited by percentiles 5 and 95 (Figure 7a), indicating that the moisture content of organic layers was not especially low during summer. Indeed, the summer of 2003 was preceded by a wet winter [16] that favored the accumulation of biomass. In turn, the time series of daily $DSR$ presents a sequence of peaks that are much beyond percentile 95, the two peaks in August being superimposed on a sequence of high values of $DSR$ throughout the month (Figure 7b). This is consistent with the extremely warm and dry weather that affected Portugal in the first two weeks of August, which was further aggravated by the very strong heatwave that struck in the first four days of that month [16]. These exceptionally warm and dry weather conditions played the double role of stressing the accumulated biomass and favoring the onset and spread of a relatively small number of very large wildfire episodes [3] that concentrated both in space (over two relatively confined regions in central and south-western Portugal) and time (within a two-week period). As opposed to the other two extreme years (2005 and 2017), 2003, the second topmost year, was not associated with prolonged drought conditions, and therefore the prevailing meteorological conditions were not adequately characterized by $CDSR$. However, it is worth stressing that such exceptional atmospheric conditions only occur once in 44 years.

The Rayleigh model with $CDSR$ as a covariate describes both the interannual variability of $BA$ and the long-term changes associated to the positive trend of $CDSR$. For instance, when excluding 2003, the interannual variability of the modeled mean explains almost 60% of the variance of $BA$. The unexplained fraction is to be attributed to other factors, namely those related to landscape management and human activities. In this respect, it is worth emphasizing the low values (below the first quartile) of the cdf of $BA$ after 2017 (Figure 6), which are likely to reflect the positive impact of policies implemented after that tragic year, aiming at better landscape management, more efficient use of prescribed burning, more organized firefighting, and a more effective reduction of man-made ignitions.

On the other hand, the positive trend in the scale parameter $\sigma$ associated with the positive trend in $CDSR$ implies long-term increases in both the mean and the standard deviation of $BA$ that, in turn, lead to a decrease in periods of return for given thresholds of $BA$. Such decreases are estimated by computing the cdf curves of $BA$ associated with the first and the last 30-year subperiods of 1980–2023. A large decrease is observed in return periods, with the relative decrease becoming larger with an increasing threshold of $BA$, from a 35% decrease for amounts greater than 120 thousand hectares up to a 59% decrease for amounts greater than 200 thousand hectares. This result emphasizes the crucial role played by climate change in the fire regimes in Portugal and severely constrains future actions regarding fire management. In this regard, statistical models as the one proposed in this work can be operationally applied to produce outlooks of annual burned area in Portugal by using $CDSR$ as derived from meteorological information provided by seasonal forecasts that are currently disseminated by a variety of weather services like the European Centre for Medium-Range Weather Forecasts (ECMWF), the US National Centers for Environmental Protection (NCEP), and national centers such as the Met Office in the United Kingdom, Météo-France in France, the Deutscher Wetterdienst (DWD) in Germany, and the Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCCl) in Italy. Along the same line, the statistical model of $BA$ can also be used to generate synthetic time series of $BA$ in future scenarios of climate based on information provided by climate models.

5. Concluding Remarks

As a major disturbance with profound impacts on the dynamics of the Earth System [65,66], vegetation fires have been the object of comprehensive studies in the last three decades [67,68,69,70,71]. Vegetation fires present spatiotemporal patterns of activity [72,73] that result from the synergy among the atmosphere that delivers the oxygen, vegetation that provides the fuel, and ignitions that ensure the initial energy required for the combustion to take place. The human factor is present in the three strands of the synergic process: ignitions are anthropogenic for the most part, and reasons are as varied as preparing fields for pastures or agriculture, hunting practices, and even negligence and arson [74,75]; landcover and land use are vastly related to socio-economic activities [76]; and climate is mostly forced by socio-economic activities and anthropogenic emissions [77].

Global patterns of vegetation fire activity are modulated by changes in temperature and precipitation related to climate variability, namely its leading modes [78,79,80]. For instance, the El-Niño–Southern Oscillation (ENSO) strongly affects fire activity in a variety of regions and ecosystems, including critical ecosystems like the tropical forests [81], whereas the North Atlantic Oscillation (NAO) has deep impacts on European vegetation dynamics [23]. Other relevant climate patterns include blocking episodes and stationary ridges, whose persistent character induces changes in temperature and precipitation that strongly impact the vegetation condition in extratropical regions [82,83,84], and equatorial Kelvin waves that are associated with tropical convection and precipitation variability as well as with the onset of the monsoon season over Africa and India [85].

Climate change has profound impacts on wildfire activity not only because of the associated global warming and changes in the water cycle but also because of the forcing of the atmospheric circulation [86,87] that leads to increases in magnitude, frequency, and extension of extreme events that favor the occurrence of large fires [88]. Recent examples are the intense droughts and associated large wildfires that have been striking the Amazonia, the Cerrado, and the Pantanal regions of Brazil [64,89,90,91]. Mediterranean Europe is also worth pointing out, not only because of its proneness to be affected by severe vegetation fires [92,93] but also because the region was identified as a primary climate change hot-spot [94,95].

Besides being a suitable descriptor of the interplay among climate, landscape, and human activities [96], wildfire activity is a useful indicator of the climatological and ecological impacts of vegetation fires, a vital source of information when setting up fire-smart forest management and fire prevention approaches [97,98,99]. The rising concern about whether global warming is accelerating [100] reinforces the need for refined estimates of wildfire activity on greenhouse gas emissions [101]. Tools like the one proposed in this study that allow quantifying the impact of climate variability and climate change on wildfire activity are therefore expected to play an increasing role.

It is, however, worth noting that the usage of statistical models of burned areas such as the one proposed in this work presupposes that the relationships among weather, vegetation, and ignitions do not change over time. This shortcoming, which is especially relevant when addressing the impact of climate change, can be circumvented in part by exploring the roles played in wildfire activity by the above-mentioned three factors for different scenarios of future climate.