Relating Spatiotemporal Patterns of Forest Fires Burned Area and Duration to Diurnal Land Surface Temperature Anomalies

Forest fires are a major source of ecosystem disturbance. Vegetation reacts to meteorological factors contributing to fire danger by reducing stomatal conductance, thus leading to an increase of canopy temperature. The latter can be detected by remote sensing measurements in the thermal infrared as a deviation of observed land surface temperature (LST) from climatological values, that is as an LST anomaly. A relationship is thus expected between LST anomalies and forest fires burned area and duration. These two characteristics are indeed controlled by a large variety of both static and dynamic factors related to topography, land cover, climate, weather (including those affecting LST) and anthropic activity. To investigate the predicting capability of remote sensing measurements, rather than constructing a comprehensive model, it would be relevant to determine whether anomalies of LST affect the probability distributions of burned area and fire duration. This research approached the outlined knowledge gap through the analysis of a dataset of forest fires in Campania (Italy) covering years 2003–2011 against estimates of LST anomaly. An LST climatology was first computed from time series of daily Aqua-MODIS LST data (product MYD11A1, collection 6) over the longest available sequence of complete annual datasets (2003–2017), through the Harmonic Analysis of Time Series (HANTS) algorithm. HANTS was also used to create individual annual models of LST data, to minimize the effect of varying observation geometry and cloud contamination on LST estimates while retaining its seasonal variation. LST anomalies where thus quantified as the difference between LST annual models and LST climatology. Fire data were intersected with LST anomaly maps to associate each fire with the LST anomaly value observed at its position on the day previous to the event. Further to this step, the closest probability distribution function describing burned area and fire duration were identified against a selection of parametric models through the maximization of the Anderson-Darling goodness-of-fit. Parameters of the identified distributions conditional to LST anomaly where then determined along their confidence intervals. Results show that in the study area log-transformed burned area is described by a normal distribution, whereas log-transformed fire duration is closer to a generalized extreme value (GEV) distribution. The parameters of these distributions conditional to LST anomaly show clear trends with increasing LST anomaly; significance of this observation was verified through a likelihood ratio test. This confirmed that LST anomaly is a covariate of both burned area and fire duration. As a consequence, it was observed that conditional probabilities of extreme events appear to increase with increasing positive deviations of LST from its climatology values. This confirms the stated hypothesis that LST anomalies affect forest fires burned area and duration and highlights the informative content of time series of LST with respect to fire danger.


Introduction
Forest fires are a source of significant ecosystem damage at global scale, as they affect the biogeochemical cycle, are a source of atmospheric emissions, alter the net carbon balance, disturb forest structure and cause long-term changes in soil properties [1][2][3][4][5]. Fires also condition anthropic activities as they threaten human lives, have a negative effect on quality of life and cause economic losses [6][7][8]. Increasing concern derives from the observation that climate change is negatively affecting spatial and temporal patterns of fire disturbance [9][10][11][12].
In Mediterranean ecosystems, prolonged droughts and heat waves create the preconditions for increases in frequency and intensity of forest fires [13,14], the underlying mechanism being the reduction of live and dead fuels moisture content as a response of the soil-plant system to increased vapor-pressure deficit [11,15,16]. Vegetation response varies with species as well as with forest structure and soil/terrain characteristics and it is determined by evapotranspiration [17][18][19]. The moisture of dead fuels, which include the organic elements of forest litter such as senescent grasses, dry leaves, small twigs and compacted organic material in the topsoil, is affected by weather variations as well and it is regulated through evaporation [20][21][22]. The moisture content of both alive and dead fuels are thus affected by weather forcing and indeed vegetation stress status has been found to be related to some meteorological drought indexes, which in turn are related to moisture content of the largest size classes of dead fuels [14,[23][24][25][26][27].
The vegetation transpiration regulation mechanism reacts to water stress conditions by reducing stomatal conductance, thus leading to an increase of canopy temperature [28][29][30]. This phenomenon can be detected by satellite measurements in the thermal infrared and has been widely used in the development of methodologies based on the satellite retrieval of land surface temperature (LST) to map vegetation stress conditions [31][32][33][34][35]. As moisture content has a direct relationship with live fuels ignitability and flames propagation [36][37][38][39], a relationship between LST and forest fires may be expected [40][41][42]. Indeed, several approaches use LST in association with optical spectral indexes of vegetation greenness or moisture content to construct physically based or empirical fire danger rating systems [43][44][45][46][47][48]. Some researchers used LST to model energy budgets [40,49,50] or to estimate heat energy of pre-ignition [51] and predict fire occurrence.
Little research was conducted to relate LST to fire characteristics. Post-fire LST was used to quantify burnt severity either alone [52][53][54] or in conjunction with optical data [55,56]. Pre-event LST was used, along with other factors, to model burned area but results were ambiguous [57]. While evidence supports the hypothesis that higher surface temperature is associated with an increased fire occurrence [58][59][60], no such relationship was previously investigated against burned area or fire duration.
Burned area is indeed controlled by a large variety of both static and dynamic factors, essentially falling into five groups: topography, such as elevation, slope, south-westness (in the northern hemisphere) or north-westness (in the southern hemisphere); land cover, including vegetation type, composition, connectivity, fuel load, pyrodiversity; climate, for example annual average daily maximum and minimum temperature; weather (including active drivers of fuel moisture) such as cumulative antecedent precipitation, wind speed, relative humidity; anthropic activity, including land development, road density, distance to settlements, fire prevention and contrasting strategies [61][62][63][64][65][66][67]. While less studied, fire duration appears to be related to similar factors [68][69][70][71]. Among these factors, only those affecting vegetation moisture are actually related to LST. To investigate the predicting capability of remote sensing measurements, rather than constructing a comprehensive model, it would Remote Sens. 2018, 10,1777 3 of 20 be relevant to determine whether an increase in LST affects the probability distributions of burned area and fire duration. Since an increase of LST would be evaluated against an LST climatology, thus implicitly implying the evaluation of a delta or anomaly, the objective would in other terms be to assess if such an anomaly is a covariate of the two named fire characteristics.
The research described in this document follows from this line of reasoning and tried to verify if a relationship linking anomalies of LST to burned area and fire duration exists, thus exploring the identified knowledge gap. The analyses were performed on the study area of Campania, Italy, for which a dataset comprising more than 8800 fire events recorded between 2003 and 2011 was made available by local authorities.
The accomplishment of the stated objective first required the definition of a method for the quantification of LST anomalies. Crucial to this step was the prior modelling of an LST climatology. Indeed, multitemporal analysis was suggested as a mean to determine seasonal minima against which to assess LST values triggering fire occurrence [58,[72][73][74]. To this purpose, the longest available time series of daily Aqua-MODIS LST data was processed with the Harmonic Analysis of Time Series (HANTS) algorithm [75,76] to construct a daily pixel-wise climatology of LST. HANTS was also used to process annual time series of LST and create cloud-and noise-free annual models of daily LST. LST anomalies were finally evaluated as the difference between the LST annual models and the LST climatology.
Further steps required the analysis of the fire data towards the identification of the closest fitting probability density function describing burned area and fire duration. The fire database was then intersected with daily maps of LST anomalies and each fire was associated with the corresponding LST anomaly value occurring at the same location on the day previous to the event. Parameters of the identified distributions conditional to LST anomaly values were determined along their confidence intervals and trends were identified [77]. Finally, probability of extreme events conditional to LST anomaly were evaluated.

Study Area
The research focused on Campania, Italy (40 • 83 N, 14 • 13 E, 13,595 km 2 , Figure 1). The interest in this region is given by its position in the middle of the Mediterranean and by the diversity of its landscape. The land cover is characterized by agricultural areas prevailing in plains and semi-natural forest-dominated areas covering hills and mountains [78]. Campania is among the most densely populated regions of Mediterranean Europe, almost all fires are triggered by human activities and fire incidence is considered high [79,80].
Remote Sens. 2018, 10, x FOR PEER REVIEW  3 of 20 an LST climatology, thus implicitly implying the evaluation of a delta or anomaly, the objective would in other terms be to assess if such an anomaly is a covariate of the two named fire characteristics.
The research described in this document follows from this line of reasoning and tried to verify if a relationship linking anomalies of LST to burned area and fire duration exists, thus exploring the identified knowledge gap. The analyses were performed on the study area of Campania, Italy, for which a dataset comprising more than 8800 fire events recorded between 2003 and 2011 was made available by local authorities.
The accomplishment of the stated objective first required the definition of a method for the quantification of LST anomalies. Crucial to this step was the prior modelling of an LST climatology. Indeed, multitemporal analysis was suggested as a mean to determine seasonal minima against which to assess LST values triggering fire occurrence [58,[72][73][74]. To this purpose, the longest available time series of daily Aqua-MODIS LST data was processed with the Harmonic Analysis of Time Series (HANTS) algorithm [75,76] to construct a daily pixel-wise climatology of LST. HANTS was also used to process annual time series of LST and create cloud-and noise-free annual models of daily LST. LST anomalies were finally evaluated as the difference between the LST annual models and the LST climatology.
Further steps required the analysis of the fire data towards the identification of the closest fitting probability density function describing burned area and fire duration. The fire database was then intersected with daily maps of LST anomalies and each fire was associated with the corresponding LST anomaly value occurring at the same location on the day previous to the event. Parameters of the identified distributions conditional to LST anomaly values were determined along their confidence intervals and trends were identified [77]. Finally, probability of extreme events conditional to LST anomaly were evaluated.

Study Area
The research focused on Campania, Italy (40°83′N, 14°13′E, 13,595 km 2 , Figure 1). The interest in this region is given by its position in the middle of the Mediterranean and by the diversity of its landscape. The land cover is characterized by agricultural areas prevailing in plains and semi-natural forest-dominated areas covering hills and mountains [78]. Campania is among the most densely populated regions of Mediterranean Europe, almost all fires are triggered by human activities and fire incidence is considered high [79,80].

MODIS LST Data
A dataset of daily gridded Aqua-MODIS LST data (product MYD11A1, collection 6) from 2003 to 2017 retrieved from the Land Processes Distributed Active Archive Centre (LP DAAC, https:// e4ftl01.cr.usgs.gov/) was used for this research. MYD11 products are generated by an angle-dependent split-window algorithm exploiting the differential atmospheric absorption in MODIS bands 31 (11 µm) and 32 (12 µm) to determine LST values from radiance measurements of clear-sky pixels. The achieved mean LST error is typically within ±0.6 K and the standard deviation of validation errors is typically less than 0.5 K [81].
Product MYD11A1 contains both diurnal (13:30) and nocturnal (1:30) LST measurements, along with corresponding quality assurance information. Preference was given to diurnal rather than nocturnal data and to Aqua-MODIS rather than Terra-MODIS, to capture canopy temperature variations due to water stress occurring at the hour of the day when maximum air temperature is approximatively achieved. Retrieved LST estimates were further masked against MYD11A1 pixel-wise quality assurance (QA) metadata and only data marked as good quality (QA bits 1,0 = 00), that is retrieved at nominal radiometric and clear-sky conditions, were retained [82,83]. However, this approach does not ensure that all cloud contaminated pixels are excluded from further processing [84].

Fire Data
A dataset of about 8800 fires officially recorded in Campania between 2003 and 2011 was provided by the Natural Resources Unit of Carabinieri, a law enforcement agency in charge of forest fires prevention, firefighting, arson investigations and prosecution and burned area inventorying. The database details the presumed date and time of fire ignition, recorded date and time of fire extinction, geographic coordinates of burned area centroid, total burned area and presumed causes. While Carabinieri actually record burnt scar perimeters on a fire-by-fire basis, according to conventional practices of field surveying with GPS receivers and desk digital cartography, these were not provided for this research. However, for the purpose of this study, this is not a source of concern on the positional accuracy of the provided centroids, as only 53 fires (0.87% of the fires in the dataset) are larger than 1 km 2 , that is of a MODIS pixel in the thermal bands.
The dataset covered a range of fire seasons that were considered safe to critical in both number of fires and total burnt area. Most fires (84%) occur in the summer season, that is June to September. About 99.8% of fires are of human origin (negligent or arson). On average, 980 fires are recorded each year, leading to the loss of more than 6160 hectares of natural areas, including 4190 ha of forests.
Fires in the database were overlaid on the CORINE Land Cover (CLC) map [85] to select those occurred in vegetated and recently burned areas only. CLC maps are produced at a nominal scale of 1:100.000, with a minimum mapping unit of 25 ha and minimum width of linear elements of 100 m and are updated every six years. Fires occurred between 2003 and 2008 were intersected with CLC 2006, while fires between 2009 and 2011 were intersected with CLC 2012. A total of more than 6100 events occurred in land cover classes reported in Table 1 were used in this research, thus excluding events recorded on agricultural land.
Observed burned area encompasses five orders of magnitude, while its average is 7.1 ha. The 95th percentile of burned area is 27.8 ha; this quantity was used as a reference for extreme events in the region. Analogously, mean fire duration is 9.4 h and the 95th percentile is 27.5 h.

Modelling Temporal Patterns of LST
The modelling of the LST climatology and of LST annual models was performed by means of the Harmonic Analysis of Time Series (HANTS) algorithm [75,76]. This method was initially proposed to fill in missing or cloudy observations and to remove outliers in time series of NDVI data by exploiting its periodic behaviour [86,87] and later used also with LST series, for example [83,88,89].
A temporal sequence of N images I(x, y, t i ), i = 1, 2, . . . , N, can be described as a Fourier series: where I(x, y, t i ) is the LST retrieved from MODIS measurements at pixel longitude x, pixel latitude y, day t i when the ith image was taken, ω j is the frequency of the jth harmonic term in the Fourier series, M is the number of frequencies of the Fourier series, a j (x, y) and ϕ j (x, y) are the amplitude and the phase of the jth harmonic term. The harmonic frequencies are integer multiples of the base frequency, that is, where L is the length of the base period. Because the zero frequency has no phase, the amplitude related to the zero frequency a 0 (x, y) is equal to the average of all N observations of I(x, y, t i ) [89]. HANTS handles the Fourier analysis as a least squares curve fitting problem within an iterative approach. In the first step, the least squares curve fitting is performed using all valid data in the series. In the second step, observations that deviate from the curve determined in the first iteration more than a pre-defined threshold (the fit error tolerance, FET) in the specified direction of rejection (lower values or higher values) are removed and the remaining data are used to compute the least square curve fitting again. The iterations are repeated until either all the remaining observations are within the FET or the number of remaining data points becomes less than the specified degree of over determinedness (DOD) [90].

Evaluation of Land Surface Temperature Anomaly
HANTS was used to decompose the time series of MODIS LST retrievals into their descriptive significant periodic components. Series comprising the first three harmonics were fit to the data with two different methods: • HANTS was executed on yearly sequences 2003-2011 of daily LST data individually to construct annual models of daily LST [82]. The objective of this approach was the removal of LST variability due to undetected cloud contamination and varying observation geometry while modelling LST annual variation. The result was a collection of new annual series of daily LST maps, one for each year being considered, computed from the identified harmonic components. These were used as representative of actual measurements.

•
The algorithm was executed on the whole 2003-2017 data set, with a base period of one year, to construct a pixel-wise daily climatology of LST [88]. The need of using this climatology as a basis for the calculation of thermal anomalies suggested its evaluation from the longest available sequence of complete annual datasets of daily MODIS LST data. The output of this process is a new series of daily LST maps computed from the identified harmonic components, representative of daily climatological values of LST.
A synthesis of the HANTS parameters adopted in the two approaches is reported in Table 2. In both, the base period is one year and the number of frequencies is set to three. The direction of outliers rejection was set as "Lo", thus leading to the removal of all data points that are more than FET lower than the fitted harmonics, according to the fact that cloud contamination in pixels causes an underestimation of LST. The development of the LST climatology is based on a more relaxed FET value, as opposed to annual models, to compensate for its inter-annual variability. The degree of over determinedness was dynamically adjusted on a per-pixel basis as the half of LST estimates marked as good quality in the QA of MYD11A1 product. Table 2. Parameters used in HANTS algorithm to pre-process LST data.

HANTS Parameters LST Annual Models LST Climatology
Length In this study, both fire size and fire duration datasets were tested against the following parametric models: A land surface temperature anomaly is hereby defined as the difference between LST annual models and the LST climatology. At a given day t i the anomaly is positive when the LST annual model is higher that the LST climatology. In this sense, LST anomalies quantify the deviation of LST from the climatology value expected in that pixel ( Figure 2). The approach of using the LST annual model rather than the actual measurements of LST hinders the detection of LST variations occurring over a short period of time. Nevertheless, it quantifies the build-up of the LST anomaly throughout the dry season while filtering the variability in LST estimates due to observation geometry, residual cloud contamination and retrieval accuracy. This in turns allows the appreciation of LST anomalies with values below the accuracies provided by the MYD11A1 algorithm and provides measurements in dates when complete cloud cover is present [88].
Daily maps of LST anomaly were produced for the entire observation period between 2003 and 2011. These maps were intersected with fire data, so that each fire was associated with the value of the LST anomaly in the same location on the day previous to the event.

Parametric Distributions of Burned Area and Fire Duration
Burned area and fire duration relate to the difficulty of control of fires and to the damage they cause and thus to fire danger [91,92]. Prior to any investigation on the relationship between these fire characteristics and LST anomalies, their probability distributions were identified.
Several parametric distributions are reported in literature to fit burned area, including normal [93], log-normal [94,95], exponential [93,96,97], gamma [97], generalized extreme value [98] and Weibull [99]. A limited number of papers report on parametric distributions of fire duration for example [100]. The diversity of these results highlights that no one single model can be identified to describe fire size distribution globally, due to the diversity of encompassed terrain, climate, ecology, forest and fire management practices [99,101] and that the closest fitting model needs to be identified on a regional basis.
In this study, both fire size and fire duration datasets were tested against the following parametric models: • Weibull These models were fitted to burned area and fire duration data by minimizing the Anderson-Darling distance [102]: where F(x) is the model cumulative distribution function and F n (x) is the empirical cumulative distribution of the sample. The maximum goodness-of-fit criterion with Anderson-Darling distance gives more weight to the tails of the distribution. The closest fitting model for each of the two variables was identified as the one providing the minimum Anderson-Darlin distance.

Conditional Distribution of Fire Characteristics
The parameters of the closest fitting distributions identified for burned area and fire duration conditional to LST anomaly were evaluated by dividing the values attained by this covariate at fire locations the day previous to the event into ten bins each corresponding to a decile, following the approach proposed in Reference [77]. In accordance with the analyses performed under the previous section, conditional parameters were determined with the maximum goodness-of-fit criterion, while their corresponding 95% confidence intervals were determined by means of 1000 bootstrap parameter estimations. Significance of the variation of observed distribution parameters across the decile bins of LST anomaly were finally verified through a likelihood ratio test where the likelihood of the model describing the entire dataset was compared against the sum of the likelihoods of the ten models in each bin.

Conditional Distribution of Fire Characteristics
The parameters of the closest fitting distributions identified for burned area and fire duration conditional to LST anomaly were evaluated by dividing the values attained by this covariate at fire locations the day previous to the event into ten bins each corresponding to a decile, following the approach proposed in Reference [77]. In accordance with the analyses performed under the previous section, conditional parameters were determined with the maximum goodness-of-fit criterion, while their corresponding 95% confidence intervals were determined by means of 1000 bootstrap parameter estimations. Significance of the variation of observed distribution parameters across the decile bins of LST anomaly were finally verified through a likelihood ratio test where the likelihood of the model describing the entire dataset was compared against the sum of the likelihoods of the ten models in each bin.

Evaluation of Land Surface Temperature Anomaly
LST in the study area exhibits significant inter-annual variability, as demonstrated by the two sample maps extracted from LST annual models of years 2007 and 2011 on the same date ( Figure 3). However, absolute LST values are not indicative of deviations from a climatology. Indeed maps of LST anomaly on the same dates show significantly different spatial patterns as opposed to the corresponding maps of LST ( Figure 4).
In the proposed sample dates, Figure 4 also reports fires occurred in the following day, represented with circles proportional to burned area. Fires occurring where a higher LST anomaly is reported result in a larger burned area, albeit such qualitative observations may vary with the chosen dates. Similar observations could be drawn for fire duration, but for the sake of brevity these are not reported herein.
Maps of LST anomaly were sampled at each fire location on the day previous to the event. Observed average LST anomaly is 1.3 K. 77% of fires occur when a positive LST anomaly is observed. On a monthly basis, this percentage varies between 69% and 88%, the only exception being December with 57%. A partial dependence of LST anomaly values from CLC classes is observed, with coniferous forest and sclerophyllous vegetation showing a wider proportion of fires occurring when a negative LST anomaly is recorded ( Figure 5).

Statistical Models of Burned Area and Fire Duration
Burned area and fire duration range over several orders of magnitude and are strongly positively skewed. For this reason, they were preliminary scaled and converted to their base 10 logarithm, so as to have log-transformed positive values only. An initial investigation highlighted that log-transformed burned area and log-transformed fire duration show a linear correlation with a In the proposed sample dates, Figure 4 also reports fires occurred in the following day, represented with circles proportional to burned area. Fires occurring where a higher LST anomaly is reported result in a larger burned area, albeit such qualitative observations may vary with the chosen dates. Similar observations could be drawn for fire duration, but for the sake of brevity these are not reported herein. Maps of LST anomaly were sampled at each fire location on the day previous to the event. Observed average LST anomaly is 1.3 K. 77% of fires occur when a positive LST anomaly is observed. On a monthly basis, this percentage varies between 69% and 88%, the only exception being December with 57%. A partial dependence of LST anomaly values from CLC classes is observed, with coniferous forest and sclerophyllous vegetation showing a wider proportion of fires occurring when a negative LST anomaly is recorded ( Figure 5). against corresponding LST anomaly. Pearson's correlation coefficients are 0.15 and 0.16 respectively, confirming that the observed large variability in these fire characteristics cannot be explained by the sole LST anomaly and indeed no trends can be clearly identified. To facilitate interpretation, data were subdivided in ten decile bins of LST anomaly and mean log-transformed burned area and logtransformed fire duration were calculated in each bin. Results plotted in Figure 8 show clearer trends, with mean log-transformed burned area increasing with increasing LST anomaly, confirming that the latter might be a covariate of this fire characteristic. A similar consideration may be drawn for logtransformed fire duration.
The parametric probability distributions identified in section 2.4. were fitted to log-transformed burned area and log-transformed fire duration using the maximum goodness-of-fit method. The corresponding Anderson-Darling distance reported in Table 3 show that log-transformed burned area is more closely described by a normal distribution, whereas log-transformed fire duration is closer to a GEV distribution. The corresponding Q-Q plots are reported in Figure 9.   against corresponding LST anomaly. Pearson's correlation coefficients are 0.15 and 0.16 respectively, confirming that the observed large variability in these fire characteristics cannot be explained by the sole LST anomaly and indeed no trends can be clearly identified. To facilitate interpretation, data were subdivided in ten decile bins of LST anomaly and mean log-transformed burned area and logtransformed fire duration were calculated in each bin. Results plotted in Figure 8 show clearer trends, with mean log-transformed burned area increasing with increasing LST anomaly, confirming that the latter might be a covariate of this fire characteristic. A similar consideration may be drawn for logtransformed fire duration.
The parametric probability distributions identified in section 2.4. were fitted to log-transformed burned area and log-transformed fire duration using the maximum goodness-of-fit method. The corresponding Anderson-Darling distance reported in Table 3 show that log-transformed burned area is more closely described by a normal distribution, whereas log-transformed fire duration is closer to a GEV distribution. The corresponding Q-Q plots are reported in Figure 9.   (Table 1). Figure 5. Boxplot of observed LST anomaly in each CLC class (Table 1).

Statistical Models of Burned Area and Fire Duration
Burned area and fire duration range over several orders of magnitude and are strongly positively skewed. For this reason, they were preliminary scaled and converted to their base 10 logarithm, so as to have log-transformed positive values only. An initial investigation highlighted that log-transformed burned area and log-transformed fire duration show a linear correlation with a Pearson's correlation coefficient of 0.57 ( Figure 6). While clearly related, the relative weakness of this correlation confirms that the two quantities are not redundant. Figure 7 reports scatterplots of log-transformed burned area and log-transformed fire duration against corresponding LST anomaly. Pearson's correlation coefficients are 0.15 and 0.16 respectively, confirming that the observed large variability in these fire characteristics cannot be explained by the sole LST anomaly and indeed no trends can be clearly identified. To facilitate interpretation, data were subdivided in ten decile bins of LST anomaly and mean log-transformed burned area and log-transformed fire duration were calculated in each bin. Results plotted in Figure 8 show clearer trends, with mean log-transformed burned area increasing with increasing LST anomaly, confirming that the latter might be a covariate of this fire characteristic. A similar consideration may be drawn for log-transformed fire duration.             The parametric probability distributions identified in Section 2.4 were fitted to log-transformed burned area and log-transformed fire duration using the maximum goodness-of-fit method. The corresponding Anderson-Darling distance reported in Table 3 show that log-transformed burned area is more closely described by a normal distribution, whereas log-transformed fire duration is closer to a GEV distribution. The corresponding Q-Q plots are reported in Figure 9.

Conditional Distribution of Burned Area and Fire Duration
Parameters of the normal distribution of log-transformed burned area show a clear trend against LST anomaly ( Figure 10). The sum of the likelihoods of the ten models fitted to burned area data in each bin was compared against the likelihood of the model describing the entire dataset by means of a likelihood ratio test. The null hypothesis in which the ten distributions are identical to the distribution describing all burned area data collectively was rejected at a significance level of p < 0.05. A similar result was observed with the GEV distribution of log-transformed fire duration, where the location a, scale b and shape s exhibit a clear trend with LST anomaly (Figure 11) and the null hypothesis is rejected with p < 0.05.
The retrieved conditional distributions (ten for each of the two fire characteristics, one in each of LST anomaly decile bins) were used to calculate the probability of fires larger than 27.8 ha (95th percentile of burned area) and the probability of fires lasting more than 27.5 hours (95th percentile of fire duration). Resulting plots ( Figure 12) show that probability of large fires ranges from 1.8% in the first LST anomaly decile to 9.9% in the tenth decile, that is when LST anomaly increases from -1.6 to 3.9 K. Analogously, probability of fires lasting more than 27.5 h ranges from 0.4% to 8.9%.

Conditional Distribution of Burned Area and Fire Duration
Parameters of the normal distribution of log-transformed burned area show a clear trend against LST anomaly ( Figure 10). The sum of the likelihoods of the ten models fitted to burned area data in each bin was compared against the likelihood of the model describing the entire dataset by means of a likelihood ratio test. The null hypothesis in which the ten distributions are identical to the distribution describing all burned area data collectively was rejected at a significance level of p < 0.05. A similar result was observed with the GEV distribution of log-transformed fire duration, where the location a, scale b and shape s exhibit a clear trend with LST anomaly (Figure 11) and the null hypothesis is rejected with p < 0.05.
The retrieved conditional distributions (ten for each of the two fire characteristics, one in each of LST anomaly decile bins) were used to calculate the probability of fires larger than 27.8 ha (95th percentile of burned area) and the probability of fires lasting more than 27.5 hours (95th percentile of fire duration). Resulting plots (Figure 12) show that probability of large fires ranges from 1.8% in the first LST anomaly decile to 9.9% in the tenth decile, that is when LST anomaly increases from -1.6 to 3.9 K. Analogously, probability of fires lasting more than 27.5 h ranges from 0.4% to 8.9%.   Figure 11. Plot of location a (a), scale b (b) and shape s (c) parameters of generalized extreme value distribution of log-transformed fire duration and their 95% confidence intervals, conditional to LST anomaly in 10 decile bins.
The retrieved conditional distributions (ten for each of the two fire characteristics, one in each of LST anomaly decile bins) were used to calculate the probability of fires larger than 27.8 ha (95th percentile of burned area) and the probability of fires lasting more than 27.5 hours (95th percentile of fire duration). Resulting plots ( Figure 12) show that probability of large fires ranges from 1.8% in the first LST anomaly decile to 9.9% in the tenth decile, that is when LST anomaly increases from −1.6 to 3.9 K. Analogously, probability of fires lasting more than 27.5 h ranges from 0.4% to 8.9%. Figure 11. Plot of location a (a), scale b (b) and shape s (c) parameters of generalized extreme value distribution of log-transformed fire duration and their 95% confidence intervals, conditional to LST anomaly in 10 decile bins.

Discussion
A clear relationship exists in several ecosystems between forest fire activity and meteorological forcing [11][12][13][14][15][16]. Prolonged absence of rainfall and increased air temperatures, while creating the preconditions for forest fires, push vegetation towards water stress conditions to which it responds by reducing transpiration. This leads in turn to an increase of vegetation temperature, a phenomenon that can be detected by remote sensing measurements in the thermal infrared [28][29][30][31][32][33][34][35]. Rather than modelling a direct dependence, this study hypothesized that remote observations of LST and more

Discussion
A clear relationship exists in several ecosystems between forest fire activity and meteorological forcing [11][12][13][14][15][16]. Prolonged absence of rainfall and increased air temperatures, while creating the preconditions for forest fires, push vegetation towards water stress conditions to which it responds by reducing transpiration. This leads in turn to an increase of vegetation temperature, a phenomenon that can be detected by remote sensing measurements in the thermal infrared [28][29][30][31][32][33][34][35]. Rather than modelling a direct dependence, this study hypothesized that remote observations of LST and more specifically deviations of LST values from a climatology, could be a covariate of burned area and fire duration of individual fires. Intuitively, prolonged events could eventually lead to larger burnt scars. Indeed, a correlation was found between log-transformed burned area and log-transformed fire duration (Figure 6), albeit its relative weakness supports the idea that these two quantities could be studied separately.
This research faced two substantial challenges. The first one was the identification of a suitable probability distribution model describing burned area and fire duration data in the study area. Several models are reported in cited literature [93][94][95][96][97][98][99][100] and indeed the diversity of these results highlights that no one single model can be identified to describe fire size distribution globally and that a model should be adopted on a per-study basis [99,101]. Among those tested herein, normal appears to be the closest fitting model for log-transformed burned area and GEV for log-transformed fire duration. In both circumstances, the fitting was not perfect towards the tails, as demonstrated by the relatively high Anderson-Darling distance (Table 3) and by the Q-Q plots ( Figure 9). Indeed, the final extent of a fire and its duration are contributed by a large number of factors related to topography, land cover, climate, weather and anthropic action. The complex and varied landscape in the study area, with significant variations of population density, topography, land use/land cover and land management practices across its extent, created a unique combination of factors shaping the probability distributions of fire characteristics that is not properly captured by the tested models.
The second challenge was the construction of a climatology of LST to use as a basis for the evaluation of the LST anomaly. The phenomenon to be indirectly detected as a deviation from the LST climatology is the reduction in stomatal conductance due to plant water stress. The need to acknowledge for its intra-annual variability excluded the opportunity to identify a seasonal mean on any base that is not over short periods of time. Indeed, the availability of daily diurnal Aqua-MODIS measurements over fifteen years allowed the calculation of an LST climatology on a daily basis. Rather than a daily average, the latter was the result of the modelling of the time series by means of the HANTS algorithm. This approach has the advantage to retain seasonal variability while filtering out disturbance sources such as undetected cloud contamination, which induce a bias by reducing the detected temperature and the varying observation geometry [89].
HANTS algorithm was also used to model annual LST series individually. This led to the definition of the LST anomaly as a deviation of the LST annual model from the LST climatology. While the use of the LST annual model instead of current measurements might hinder the detection of LST variations occurring over a short period of time, this approach has the advantage of still detecting the build-up of the LST anomaly throughout the dry season while leveraging the named effects of observation geometry and residual cloud contamination [88].
The LST climatology and the LST annual models were constructed using the same HANTS parameters, with the only exception of FET ( Table 2). The need for a more relaxed FET in the LST climatology was justified by the need to account for the inter-annual variability of LST. While affecting the shape of the modelled curve, a lower FET would clearly determine a climatology characterized by generally higher LST values. The value of FET = 6 K was identified as a compromise between the need of constructing an LST climatology of general validity while still rejecting cloud-contaminated data in both cooler and hotter years.
Fire data were intersected with maps of LST anomaly and each fire was associated with the LST anomaly value recorded at fire location on the day previous to the event. The underlying idea is that LST data is produced in near-real-time by ground receiving stations, allowing the mapping of fire danger forecast for the following day. The way LST anomaly was defined in this paper, that is as a deviation of LST annual model from LST climatology as opposed to the deviation of the actual MYD11A1 estimate, implies a slow day-to-day variation of LST anomaly. This in turn increases the temporal validity of produced LST anomaly maps up to a certain extent. Checking the effect of sampling in time at more days before the event was beyond the objectives of this research. However it is here anticipated that tests performed on LST anomalies recorded five days before the event led to results similar to those reported herein.
Observed average LST anomaly is 1.3 K, while 23% of fires occur with a negative LST anomaly. The distribution of LST anomaly appears to be partially dependent on land cover class as reported in Figure 5. Negative values are proportionally more prevalent in coniferous forest and sclerophyllous vegetation than in other land cover classes. The observation of this dependence was expected from literature review (e.g., [63,65]) and does not affect the quality of further findings. Indeed, evaluations of burned area and fire duration were performed conditional to LST anomaly, that is leveraging out all other parameters. In effect, other factors such as accessibility of the zone, availability and efficacy of the fire extinguishing means and the winds can heavily influence the occurrence of fires, the final burned area and event duration regardless of the previous LST anomaly, either positive or negative This study demonstrated the informative content of time series of LST through the observation that LST anomaly is a covariate of burned area and fire duration. Fire data were grouped in ten decile bins of the associated LST anomaly values and the parameters of the identified distributions along with their 95% confidence interval were evaluated in each of them. The choice of the number of bins was initially tested with a trial and error approach, towards the identification of a compromise between the clarity of observed trends and the amplitude of the distribution parameters confidence interval. As similar results were observed across a range from five to twenty bins, for the sake of objectiveness in the approach the number of ten was chosen as the most appropriate.
Mean and variance of the normal distribution describing log-transformed burned area both tend to increase with increasing LST anomaly values (Figure 10), likewise the location, scale and shape parameters of the GEV distribution describing log-transformed fire duration (Figure 11). A likelihood ratio test confirmed that probability distribution models of burned area and fire duration conditional to LST anomaly are collectively significantly different than the models describing the entire dataset. Along with the observed trends in parameters values, this result confirms the stated hypothesis that LST anomaly as defined in this paper is a meaningful variable contributing to fire danger, that is that this quantity, along with its time-dependent nature, may be used pairwise other relevant parameters towards the statistical modelling of burned area and fire duration.
As a consequence of the variation of burned area and fire duration probability distributions conditional to LST anomaly, it is possible to plot and interpret how the probability of extreme events evolves with increasing LST anomaly values. Indeed, end users such as forest managers and civil protection agencies are particularly interested in these probabilities to drive their preparedness activities [103]. The models constructed in the ten decile bins of LST anomaly demonstrate that the probability of a forest fire to result in a total burned area exceeding a given threshold (in our example, the 95th percentile of the dataset) significantly increases with increasing LST values (Figure 12a). A similar result was found for fire duration (Figure 12b), supporting the idea that maps of LST anomaly such as those in Figure 4 are useful to depict the contribution of LST to fire danger.

Conclusions
Vegetation response to meteorological factors contributing to fire danger-prolonged absence of rainfall and high air temperature-results in an increase of LST that can be detected by remote sensing measurements in the thermal infrared as a deviation from climatological values. This paper demonstrates that such LST anomalies are a covariate of forest fires burned area and duration. While several studies demonstrated how a wide number of both static and dynamic factors related to topography, land cover, climate, weather and anthropic activity affect the probability distribution of these two fire characteristics, to the best of authors' knowledge no previous research was conducted to investigate the role of satellite measurements of LST.
The initial hypothesis was addressed by first identifying probability distributions functions describing available fire data. Among those tested, log-transformed burned area is closer to a normal distribution, while log-transformed fire duration is closer to a generalized extreme value distribution. The HANTS algorithm was then used to process time series of diurnal Aqua-MODIS LST measurements and construct a climatology against which anomalies of LST were quantified. Finally, parameters of the identified distributions conditional to LST anomaly where then evaluated, showing clear trends.
The observed variation of burned area and fire duration distributions conditional to LST anomaly demonstrate that increasing positive deviations of LST from the expected seasonal value correspond to an increasing probability of extreme events, that is, of the final fire extent and duration exceeding a given threshold. This observation clearly identifies a practical mean to interpret maps of LST anomaly. As opposed to typical fire danger rating systems based on meteorological data, this remote sensing quantity has the advantage of claiming a higher spatial resolution. It is here highlighted that the identified relationships are preconditioning in nature and do not predict actual fire occurrence. The latter is related to a different array of determinants relating to probability of a heat source leading to an ignition.
This analysis was performed ex post trough the evaluation of LST annual models against an LST climatology. This approach is observational in nature and not predictive. While the achievement of more generality of the proposed research would require investigations in a wider and diverse array of regions, the predictive capability of the developed approach still needs to be demonstrated. This was beyond the objectives of this research and will be the subject of further investigations.