Analysis of the Relationship between Land Surface Temperature and Wildfire Severity in a Series of Landsat Images

The paper assesses spatio-temporal patterns of land surface temperature (LST) and fire severity in the Las Hurdes wildfire of Pinus pinaster forest, which occurred in July 2009, in Extremadura (Spain), from a time series of fifteen Landsat 5 TM images corresponding to 27 post-fire months. The differenced Normalized Burn Ratio (dNBR) was used to evaluate burn severity. The mono-window algorithm was applied to estimate LST from the Landsat thermal band. The burned zones underwent a significant increase in LST after fire. Statistically significant differences have been detected between the LST within regions of burn severity categories. More substantial changes in LST are observed in zones of greater fire severity, which can be explained by the lower emissivity of combustion products found in the burned area and changes in the energy balance related to vegetation removal. As time progresses over the 27 months after fire, LST differences decrease due to vegetation regeneration. The differences in LST and Normalized Difference Vegetation Index (NDVI) values between burn severity categories in each image are highly correlated (r = 0.84). Spatial patterns of severity and post-fire LST obtained from Landsat time series enable an evaluation of the relationship between these variables to predict the natural dynamics of burned areas.


Introduction
Land surface temperature (LST) is one of the most important factors controlling physical processes responsible for the land surface balance of water, energy and CO 2 [1][2][3].In the context of wildfire studies, fire-induced environmental changes cause variations in the spatial distribution of LST, mainly due to a decrease in transpiration and an increase in the Bowen ratio (β = sensible heating/latent heating) [4].The higher post-fire LST of the burned areas was observed in field data [5,6] and remotely-sensed images [7,8].Moreover, according to Beringer et al. [9], there is a relationship between fire intensity and an increase in the Bowen ratio, as far as fire intensity determines the likely impact on energy and carbon fluxes.Consequently, burn severity, defined for the current study as the amount of change in a burned area with respect to the pre-fire conditions [10][11][12], is very dependent on fire intensity [13] and can be considered a key variable in understanding the spatial distribution of LST in the immediate post-fire environment [8].The regrowth of vegetation is also one of the most important factors controlling LST in the years following a fire, as vegetation cover and bare ground have different emissivity, defined as the ratio between the object emitting capacity and that of a blackbody at the same temperature.That is why spatio-temporal patterns of LST can help monitor the processes that structure ecosystem development and may assist in developing appropriate management strategies following forest fires.
Satellite sensors have long been used in wildfire research [11,14] to assess variables related to burn severity and vegetation recovery in a cost-effective and time-efficient way (among others [15,16].On the medium spatial scale, Landsat has provided global coverage since 1984, with Landsat 8 launched at the beginning of 2013, ensuring the continuity of data record [17].However, compared to optical bands, the use of Landsat thermal data presents additional challenges [1,18]. Radiance levels in the thermal region of the spectrum depend not only on the amount of solar radiation received, but also on the ability of the surface to emit energy, expressed by its emissivity and atmospheric conditions (water vapor and temperature).At present, several physically-based methods have been suggested for LST estimation based on thermal infrared data from satellites, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) and Landsat [19].In the specific case of Landsat-5 TM, with only one thermal infrared band available, atmospheric profiles of temperature and water vapor content must be known for the exact time of image acquisition, as well as the knowledge of surface emissivity for each pixel.There are several LST algorithms applicable to Landsat 5/7, including mono-window [20], single channel [21,22] and the on-line Atmospheric Correction Parameters Calculator (ACPC) [23,24].All of the procedures report similar estimation errors of 1-2 K.
Burn severity can be assessed through the calculation of spectral indices, which are focused on reflectance changes in burned areas mainly related to vegetation removal, soil exposure, changes in water content and the deposition of carbon and ash [25].Although the Normalized Difference Vegetation Index (NDVI) [26] yields good results for burn severity assessment [27,28], the Normalized Burn Ratio applied in a two-date approach, the delta Normalized Burn Ratio (dNBR) [12], outperforms the other indices [29][30][31][32][33]. dNBR can be considered a consistent method for burn severity assessment, due to its proven relationship with field severity metrics.Empirical models have shown strong relationships (r 2 > 0.6-0.7) between dNBR and specific parameters of burn severity, such as ash cover percentage, tree mortality, or twig diameter [30,[34][35][36][37], or field indices, such as the Composite Burn Index (CBI) [31,[38][39][40][41][42].Moreover, the bi-temporal approach, where values of the post-fire image are subtracted from values of the pre-fire image, is considered the best approach to detect change caused by fire.Spectral vegetation indices have been proven useful in monitoring seasonal variations in vegetation development (phenological cycle) [43,44], as well as post-fire plant regeneration [45,46]: strong correlations were observed between the NDVI and various biophysical vegetation parameters, such as Leaf Area Index (LAI), the fraction of photosynthetically active radiation (fPAR) or vegetation abundance [47].
Although relationships between burn severity, NDVI and LST values seem quite clear, few studies have explored these [8,48,49].There are indications that the inclusion of thermal information in spectral indices for severity mapping improves their performance [48,49].The post-fire LST-severity relationship was assessed by Veraverbeke et al. [8] using MODIS images for a two-year period after fire, detecting an increase in post-fire LST up to 8.4 °C for a conifer forest.However, Landsat images can be especially suitable, because both the severity and LST of burned areas can be estimated in a more detailed spatial resolution.Therefore, the objectives of this study are: (1) to evaluate changes in LST for several images over a two-year period after fire; (2) to analyze the relationship between LST and burn severity estimated using the dNBR index; and (3) to study the relationship between vegetation regrowth measured by NDVI and changes in LST.The working hypothesis tested in this study is that the spatial distribution of LST in the burned areas depends on burn severity and that the LST range in each image is related to the phenological cycle and the time elapsed since the fire.From a methodological perspective, this study relies on the potential of remotely-sensed data and, more specifically, Landsat data to estimate LST, burn severity and vegetation regrowth.

Study Area
The study area of the Las Hurdes 2009 wildfire is located in Extremadura, in the province of Cá ceres, Spain (40°19ʹ-40°24ʹN, 6°10ʹ-6°15ʹW) (Figure 1).It is a hilly area with elevation ranging from 390 to 1280 m above the sea level.The typical acid fine-textured soils are mainly umbricLeptosols and humicCambisols formed over metamorphic bedrock [50].The Mediterranean climate (Csa according to the Köppen classification), characterized by an annual average temperature of 16°C and approximately 550 mm of precipitation, has a four-month hot, dry period from June to September [51].
The Las Hurdes fire analyzed in this study burned more than 3000 ha of the 30-40 year-old pine forest (Pinuspinaster) in four days (25-28 July 2009).According to the Spanish Third National Forest Inventory [50] (sample points shown as points in Figure 1), the average tree coverage is around 40%; besides Pinuspinaster, otherspecies, notablyArbutus unedoandQuercus ilex are also present.In Spain, Pinus pinaster occupies more than 1 million ha and is highly important to Spanish forestry [52].It is also the species most affected by wildfires (27.96% of the burned area) [50].Growth usually occurs in spring (early April to mid-June) and autumn (late August to early October) [53].The seed production is generally related to the fire regime.Stands suffering recurrent, high-intensity fires show more serotinous cones and a large aerial seed bank compared to stands where crown fires are not frequent [54].

Data
Data from recently calibrated Landsat-5 TM archive [55] were used in this study.Landsat-5 TM images are composed of six optical and one thermal (bandwidth of 10.4-12.5 μm) spectral bands.Spatial resolution is 30 m for optical bands and 120 m for the thermal band.
Fifteen clear sky images, path 202/row 32, covering the period from July 2009 to September 2011, downloaded from the NASA website [56], are listed in Table 1 along with the information on the observation geometry and atmospheric conditions (near-surface air temperature T air and relative humidity We used preprocessed level L1T Landsatdata.The downloaded images (GeoTiff format) were available in the UTM projection (datum: WGS84).The digital elevation model with 25-m resolution in the UTM projection was downloaded from the online archive of the National Center for Geographic Information (Spain) [59].It was processed using ArcGIS software [60] to obtain information on the surface slope and aspect.

Atmospheric Correction of the Optical Bands
An open-source Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) from NASA Goddard Space Flight Center (GSFC) [61] was used for the atmospheric correction of the optical bands.It obtains parameters required for atmospheric correction from the National Centers for Environmental Prediction (NCEP) reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, on-line [62](atmospheric pressure and water vapor), at 2.5° spatial resolution and the Earth Probe Total Ozone Mapping Spectrometer (EP TOMS) (ozone) at 1° spatial resolution, available from [63].The obtained values are resampled to the same spatial resolution of 1.2 km and each image is processed and corrected independently.One of the advantages of the system compared to other similar tools, is that it takes the original Landsat data (DN values) as inputs and provides atmospherically corrected reflectance values for each of the optical bands as outputs ,without the need for previous data transformation or scaling by the user.

Monitoring Vegetation Recovery
Monitoring of vegetation recovery was performed through NDVI calculated for each available image.The NDVI is based on the difference between the maximum reflection of radiation in the near-infrared spectral bands (0.78-0.90 μm) and the maximum absorption of radiation in the red spectral band (0.63-0.69 μm).The difference of the reflectances is normalized by their sum, reducing the effect of shadows, resulting in NDVI = (NIR − VIS)/(NIR + VIS).
Values of the NDVI range between −1.0 and +1.0.The wide use of NDVI for vegetation monitoring arises because of its positive correlation with characteristics of plant status and abundance.NDVI frequently serves as a proxy for biomass, although the relationship between them is often non-linear [26,44],and NDVI shows saturation before biomass reaches its maximum levels.In spite of the limitations, NDVI is commonly used in assessing vegetation recovery after fire (among others [8,16,27,64]).It is sometimes used as a metric of burn severity [8,48,65,66].

LST Estimation
LST was calculated using the mono-window (Mw) method [20].Prior to the LST estimation, band 6 original data were transformed first into radiance, with the help of the data from the header files, and next into the at-sensor brightness temperature.The Mw algorithm [20] requires three a priori known parameters: atmospheric transmissivity (τ) calculated from the water vapor content, effective mean atmospheric temperature (T a ) and surface emissivity (ε).The formula used to calculate LST (Ts) is the approximation of the radiative transfer formula and includes two empirical coefficients a and b: (1) where a = −67.355351and b = 0.458606 are constants, T sensor is the at-sensor brightness temperature and: Formulas for the estimation of the atmospheric correction parameters were developed by Qin and Karnieli [20] using LOWTRAN 7 simulations.The simulation of atmospheric transmissivity τ, depending on water vapor content, yielded Equation (3a,b) for a low temperature profile (18°C) and Equation (3c,d) for a high temperature profile (35°C) [20]: The effective mean temperature T a is computed for specific atmospheric conditions using Formulas 4a-c based on the ratio of water vapor content at a particular altitude to total atmospheric water vapor content and near-surface local air temperature T 0 [20]: T a = 17.9769 + 0.91715T 0 (tropical atmosphere) (4c) The empirical formula (Equation ( 5)) developed by Butler [67] based on Bolton [68] and adjusted for central Spain by De Vicente and Pulido [69] was used to estimate atmospheric water vapor content: (5) where w is the water vapor content (g• cm −2 ), T 0 is the near-surface air temperature in °C and RH is the relative humidity (%).
When working with Landsat thermal data, surface emissivity estimation required for calculating LST is a challenge, because only one thermal band is available.To solve the problem, the NDVI-based methods, which rely on the information from the image used for the LST retrieval, were successfully applied [70].One of these, the NDVI thresholds method (NDVI THM ) [71,72], based on the findings of Valor and Caselles [73], was used to calculate surface emissivity in this study.The emissivity for different NDVI ranges was estimated using different functions.For water and fully vegetated pixels, the emissivity values of 0.985 and 0.99, respectively, were assigned following the suggestion of Sobrino et al. [74].The soil emissivity value of 0.984 is a result of the field measurements using the box method [75] and is similar to values reported by previous research [74,76].As for the pixels with the mixed cover of vegetation and soil (0.1 ≤ NDVI ≤ 0.7), emissivity ε is calculated using Formula (6) [72-74], which involves vegetation fraction P V estimated from a scaled NDVI, according to Choudhury et al. [77] and Gutman and Ignatov [78] (Equation ( 7)): 0.990 0.984 1 0. ( ) ) 1 ( 04 where NDVI pixel is the NDVI value of a pixel.NDVI thresholds for the mixed pixels range are based on image histogram analysis.

Burn Severity Estimation
In this research, dNBR was the spectral index applied for burn severity evaluation due to the very strong association observed between dNBR and field burn severity measurements in conifer forests [34,42,79] and more specifically, in Mediterranean areas [80,81].Likewise, as LST values were obtained from Landsat data, it was considered appropriate to use the burn severity index especially designed for Landsat spatial and spectral specifications [12].The methodology followed for dNBR calculation was [82]: (1) pre-and post-fire images were transformed to reflectance R and atmospherically corrected; (2) an NBR image was generated for both dates using the formula (R 4 − R 7 )/(R 4 + R 7 ), where subscripts correspond to the band numbers; (3) dNBR was calculated as NBRpre-fire − NBRpost-fire; and (4) the polygon encompassing fire-affected pixels (dNBR > 100) plus a 350-m buffer was defined for the purposes of analysis.
dNBR values are sometimes grouped into discrete classes of burn severity (e.g., low, moderate and high) [12].Original thresholds for these intervals were not thought to be used as fixed values, valid worldwide.Several studies have used the relationship between dNBR and CBI to calculate dNBR thresholds representing breaks between burn severity classes [39,40,83], sometimes with fairly insignificant differences from the initially suggested values [84].However, there are also studies that have adopted them in ecosystems quite different from the one for which they were created [30,31,34].For simplicity and objectivity, thedNBR values suggested by Key and Benson [12] were used for creating the burn severity categories as follows: unburned (UB) (from −100 to 99), low severity (LS) (from 100 to 269), moderate-low severity (MLS) (from 270 to 439), moderate-high severity (MHS) (from 440 to 659) and high severity (HS) (from 660 to 1300).

Statistical Procedures
The comparison of pre-and post-fire images suffers from problems related to interannual phenological differences and time since fire [12,45,85,86]: the overall regeneration trend may vary significantly from one year to another due to climatic differences.To solve this problem, Dí az-Delgado and Pons [45] compared burned and unburned plots within the same image, while Veraverbeke et al. [86] used a control plot selection procedure based on Lhermite et al. [87], which exploits the similarity between the temporal evolution of the burned and unburned pixels.In this context, two different approaches to the temporal study of the LST-severity relationship were applied in this research.First, variations in LST and NDVI differences throughout the 27 months after the fire were identified by comparing the images captured at similar moments of the annual phenological cycle in different post-fire years.This analysis was applied to images satisfying the following criteria: (1) post-fire images from different years can be compared only if the acquisition day corresponds to the same phenological stage of Pinuspinaster (all of the images used for comparison in this study are acquired within the period between two active growth phenological stages between mid-June and late August [53]);and (2) the difference in atmospheric temperature between compared dates has to be lower than 1.5°C (Table 1).Thus, the following raster arithmetic calculations were applied: (1)  Second, statistical differences were studied between the LST and NDVI values observed in the burn severity categories.To reduce the spatial auto-correlation effects, a random sample of 10% pixels by severity category, including the UB category for reference, was extracted from the pixels inside the study site perimeter (n = 4230).Sample points were analyzed independently for each date using ANOVA analysis and Tamhane'sT2 post hoc test algorithms.Moreover, for further study of the temporal differences between burn severity categories, the variables -fire severity differences in LST‖ (fsdLST) and -fire severity differences in the NDVI‖ (fsdNDVI) were analyzed.fsdLST specifically refers to the LST differences between areas within burn severity categories: UB, LS, MLS, MHS and HS.It was accomplished by Formula ( 8): (8) where is the mean value of the LST variable and i, j are a pair of burn severity categories.A similar procedure (Equation ( 9)) was applied to calculate fsdNDVI: where is the mean value of the NDVI variable and i, j are a pair of burn severity categories.

Spatial Pattern of dNBR
The RGB 7-4-3 band combination (Figure 1) depicts the Las Hurdes fire perimeter in shades of red associated with the low reflectance in the NIR band, a characteristic of zones of scarce vegetation, and high reflectance at 2.1 µm in the SWIR spectral region, typical of areas with a low moisture content.This is the typical spectral response of burned areas [79] (Figure 2).Different exposure time and different fire intensity result in the great spatial variability of burn severity in the affected ecosystem.The spatial distribution of burn severity, classified from the original dNBR threshold values, can be seen in Figure 2. Within the Las Hurdes fire, 32.9% of the burned surface presents HS, 37.4% MHS, 18% MLS and 11.7% LS.On the whole, Las Hurdes was a high severity fire, since more than 70% of the area falls within the MHS and HS categories.However, within the fire perimeter, two wide diagonals of low severity pixels divide the burned area in the north and south (Figure 2), defining four sectors: two in the north with a large number of high-severity nuclei, a very large one in the center and one of predominantly moderate-low severity in the south.The predominance of the highest burn severity intervals is also related to the initial approach applied to the burn severity assessment, by using an immediate post-fire image and not giving time for the ecosystem to show additional responses to fire [12].

Temporal Dynamics of LST and NDVI Values
This section presents the temporal dynamics of LST and NDVI throughout the study period.Descriptive statistics for LST and NDVI (Tables 2 and 3) refer to data from all of the available images: the pre-fire image (13 July 2009) and 14 post-fire images taken between July 2009 (one day after fire), and September 2011 (two years after fire), while Figure 3 shows data in the form of graphics on four different dates: 13 days before the fire on13 July 2009, and on three midsummer dates corresponding to successive post-fire summer seasons (29 July 2009,16 July 2010, and 4 August 2011).Values are grouped by severity categories.In addition, Figure 4 shows the spatial distribution of LST and NDVI on the same dates as Figure 3.
In the pre-fire image, all burn severity categories present similar average LST values (~30 °C) (Table 2).The coolest areas associated with greater biomass are those registering the highest severity levels after fire (Figures 3 and 4).The existence of this type of relationship between pre-fire biomass and further burn severity was previously reported by Garcí a-Martin et al. [88], who demonstrated that knowledge of crown biomass enables the prediction of the burn severity levels.The immediate effects of the fire on the LST are reflected in the first two post-fire images (29 July and 30 August 2009) closest to the event.For the visual assessment of these effects, Figure 5 presents the spatial distribution of the dLST, where LST values of the post-fire image are subtracted from the pre-fire image.According to this formula and the assigned colors, areas with the greatest increase of LST are highlighted in red, and areas without a change are shown in green.To improve the understanding of pre-and post-fire LST changes, Figure 6 presents the confidence levels of average values for the dLST by burn severity category.The average LST increase is 13°C, reaching 20°C for the HS pixels.The generalized LST increase in the post-fire image (both in the burned and unburned areas) may be due to the fact that this image was acquired on a date closer to the middle of summer than the pre-fire image, and therefore, the air temperature was high.However, thermal differences between HS and UB categories within the post-fire image (>10°C) reveal the influence of burn severity on the spatial distribution of LST (Table 3).The decrease of aboveground green biomass in the burned zones [12], especially in those of higher severity, and the appearance of lower emissivity coverage (ash, char and mineral soils) lead to a large increase in the LST.Elevated LST after fire events is mentioned by several authors (among others Lambin et al. [7]; Montes-Helu et al. [5]; Wendt et al. [6]).Veraverbeke et al. [8] studied this increase using MODIS imagery following the major Peloponnese fire in 2007.Until now, few studies have analyzed spatiotemporal patterns of post-fire surface temperature using Landsat data, although the high potential of existing single channel algorithms, such as the mono-window (MW) method by Qin et al. [20] or the single-channel (SC) method by Jimé nez-Muñoz and Sobrino [22], has already been demonstrated [89][90][91].The greater surface heterogeneity of the burned areas due to the incorporation of combustion products, changes to lighter-colored soil and ash, char and scorched, then blackened, vegetation [12] results in an increase in post-fire thermal variability (SD values ~5, Table 2).Some interesting ideas arise regarding the influence of burn severity on the LST distribution and the contrast between areas of different burn severity categories.The general decrease in LST observed in the first post-fire autumn (September and October data of 2009) is probably associated with lower solar illumination angles.When the sun is directly above the observation location and the sunlight is perpendicular to the land surface, the amount of solar radiation received by the surface is at its maximum.However, as the angle between the sun and a surface is continually changing, the surface gets only part of the incident sunlight.Topography (slope) and sun azimuth also affect the incidence angle of sunrays and the time the area is illuminated by the sun.However, burn severity remains the main factor influencing the spatial distribution of LST: higher LST values correspond to higher burn severity and vice versa (Table 2).Likewise, post-fire thermal variability within burn severity categories maintains the level observed in the immediate post-fire image.
The same patterns of LST changes are observed in the images from different years: (i) same season LST values (spring, summer, autumn) become lower from year to year; (ii) the spatial distribution of LST values is qualitatively in agreement with the burn severity categories; and (iii) differences between extreme severity categories in 2010 are slightly lower compared to 2009 (the year of fire) and even lower in 2011 (Table 2).This smoothing of contrast among categories can be explained by both the effects of time on the combustion products and, most of all, the effects of vegetation regeneration, reflected in NDVI values registered in all of the temporal series in the various burn severity categories (Table 3) and a visual comparison of three dNDVI images (Figure 7): each image is calculated as the subtraction of the NDVI raster of one of the post-fire summer seasons from the pre-fire NDVI (dNDVI 2009 , dNDVI 2010 and dNDVI 2011 ).The images show areas with higher dNDVI in red and those with lower dNDVI (similar to the pre-fire situation) in green.The progress of vegetation recovery is quite evident: while the highest differences are characteristic to the first post-fire summer, the contrast between burned and unburned areas is smoothed in 2010 and especially in 2011 data, with much lower dNDVI values in the corresponding images.This successful regeneration process is explained by the efficient recovery mechanisms of the vegetation species dominant in this area.Pinus pinaster, the main species affected by the Las Hurdes fire, is highly adapted to fire-prone environments through the massive release of seeds from serotinous cones after fire forgermination [92,93].In the same way, shrubland species observed in fire affected areas near the study site (Erica arborea, Erica lusitanica, Cistus ladanifer, Phillies angustifolia, Cytisus scoparius, Calluna vulgaris and Lavandulastoechas) [94]also apply efficient post-fire reproduction strategies in recolonizing burned areas.
Our results demonstrate that the LST increase in fire-affected areas was evident in the analyzed series of images, which cover all of the seasons of the two post-fire years, except winter.This is similar to the results reported by the previous research [5,6,8], although the range and the size of the differences between severity categories of the same date is much larger than that detected in the earlier studies.This is probably due to the different response of the analyzed vegetation: much more homogeneous in this study (predominantly conifer forests) than analyzed in the study by Veraverbeke et al. [8] (shrublands, olive groves, coniferous and deciduous forests).Immediate post-fire increment in LST calculated from Landsat is much more pronounced than that registered for similar vegetation cover at the same phenological stage registered by MODIS, because of the difference in spatial resolution between the two sensors, i.e., the smoothing of contrasts in the lower resolution images.

Analysis of fsdLST and fsdNDVI
The results of the date-by-date LST-NDVI comparison by severity categories are shown below.It can be seen that the differences between burned/unburned areas increase with burn severity in terms of LST from around 3 °C to almost 7 °C (Figure 8a) and from 0.09 to 0.21 in terms of the NDVI (Figure 8b).Mean fsdLST between the successive severity categories is 1.42°C.In fact, significant statistical differences (p < 0.05) were registered in all of the pairs, with the only exception being HS-MHS and MHS-MLS in the images from August and September of 2011.
A detailed date-by-date analysis of the differences between severity categories allows for the identification of common features.Each pair in Figure 8 shows the distribution of the fsdLST and fsdNDVI by date.The size of the circle reflects the between-category distance for each pair (i.e., a size of four corresponds to combinations of the extreme burn severity categories UB-HS).The color in these figures represents the type of categories paired: green when one of the categories in the pair is the UB; red and orange when the HS category is involved and blue for the combinations of the intermediate categories.Pairs combining high burn severity levels (HS and MHS) and the UB class register the most pronounced differences (between 7 °C and 10 °C).A seasonal pattern is observed during 2009 and 2010, as well as a stronger decrease and temporal stabilization in 2011, two years after the fire.Pairs formed by consecutive categories (HS-MHS and MHS-MLS) (size = 1) show lower differences (<3 °C), without any specific temporal pattern.Differences below 1°C are almost exclusively observed in combinations of high severity levels (HS-MHS) on all of the dates, except in the image taken just after the fire, where they are slightly above 1 °C.Between these two groups of high and low differences, differences for the HS-LS (orange), UB-MLS (light green) and MHS-LS (blue) (size = 2-3) pairs show a large range of variation (from 2 °C to 7 °C).
The greatest fsdNDVI between 0.25 and 0.35 (Figure 8b) are always related to the comparison between UB and all of the other severity categories (green).They are mainly observed in the images of the first post-fire summer, when the effects of fire are more obvious, and especially in March and April 2010.Lower fsdNDVI (0.25-0.15) are characteristic of the pairs formed by the UB and HS (size = 4) in the last images of 2010, MHS and MLS (light blue, size = 1) and also between HS-LS (orange, size = 3) until the spring of 2010.Many pairs register differences between 0.15 and 0.05.The majority of images included in this group are from 2011.Differences below 0.05 correspond to the HS-MHS pair on all of the dates and the HS-MLS pair on the dates after June 2010.
Analysis of the fsdLST and fsdNDVI in 2011 reveals: (1) lower fsdLST compared to 2010; and (2) progressive smoothing of contrast between severity categories.fsdNDVI are below 0.10 (Figure 8b), and fsdLST values are less than 5°C (Figure 8a), except in May, when they are slightly higher.A general downward trend is observed in both the fsdLST and fsdNDVI throughout the time series, especially significant in 2011.Thus, the scatterplot in Figure 9 highlights the strong relationship between these two variables (r = 0.84): most of the bigger circles in the graphic are located in the upper part of the scatterplot, except those corresponding to the 2011 dates (in green), which are located in the lower part of the plot, always below the reference line, due to the minimizing effects of vegetation regeneration on fsdLST.Unusually low fsdLST values observed in March of 2010 are due to particularly low air temperature on this date.However, the spatial distribution of LST depends not only on burn severity and its interaction with local-scale variables, such as surface emissivity (vegetation regeneration).Factors explaining intra-and inter-annual LST changes also include illumination geometry controlled by solar azimuth and elevation angles and topography reflected in slope and aspect.The role of aspect in the spatial distribution of LST is shown in Figure 10.The figure presents mean LST values for eightcategories of aspects (at 45-degree intervals) grouped by burn severity levels from the pre-fire image up to the image taken 27 months after the fire in September of 2011.
At first glance, some influence of aspect on the spatial distribution of LST and its relationship with severity levels, cover type and day of the year is observed.High LST values are systematically registered on SE-facing slopes (between 90° and 180°), with slight variations depending on the image date.Conversely, values corresponding to pixels in NNE-and NW-facing slopes (between 270° and 360° and between 0° and 45°, respectively) always register lower LST.Differences between hot and cold orientations deepen in the spring and autumn images, due to the lower elevation angles of the sun.For example, in the image from July2010, the differences between the aspect intervals described above can exceed 6°C, and in October2009, they can be higher than 15°C, because of the deeply shaded areas.However, thermal contrast between pixels of different aspects is not as pronounced in the pre-fire image, where it does not exceed 4°C.Therefore, the fire and the consequent vegetation removal lead to greater thermal heterogeneity, increasing the role of aspect in the spatial distribution of LST.In later images r = 0.84 (those from 2011), however, vegetation regeneration reduces the differences in LST between burn severity categories and aspect intervals, as can be appreciated in the 2011 images.

Conclusions
The study quantifies fire-induced changes in the spatial and temporal distribution of land surface temperature (LST) within the Las Hurdes fire in Central Spain using Landsat imagery.Immediately after the fire, the burned zones were, on average, 7.6 °C warmer than the unburned; the difference with the unburned areas was above 10 °C for the zones of high burn severity.The size of LST differences was directly related to the area's burn severity.After the first months following the fire, LST contrasts between burned and unburned areas in the same image decrease, although the LST differences between areas of different burn severity categories are still detectable for two years after fire.
The spatial distribution of post-fire LST is mainly influenced by vegetation regeneration.In the specific case of the Las Hurdes fire, study results point to the vegetation regrowth two years after the fire as a key factor in the temporal evolution of LST values, making less noticeable the consequences of fire as time elapses.
LST contrasts in the areas of different burn severity are also enhanced by the aspect and illumination geometry, being higher for the better-illuminated slopes.As vegetation recovers, the differences between aspect intervals considerably decrease.
The study draws upon the Landsat potential to provide spatially continuous quantitative estimation of land surface temperature and demonstrates the influence of burn severity and post-fire vegetation recovery, both of them assessed by spectral indices, on the spatial distribution of land surface temperature, one of the key parameters controlling physical processes in fire-altered areas.
Future research will approach the relationships between burn severity, LST and vegetation regeneration in other ecosystems and test the possibility of combining LST with commonly used metrics to improve burn severity differentiation.

Figure 1 .
Figure 1.Map of the fire site.Points indicate the location of the Spanish National Forest Inventory parcels [50].

Figure 3 .
Figure 3.The relationship between burn severity categories and LST in °C (left panel), and NDVI (right panel).Bars indicate confidence interval of average values (α = 0.01).Each graphic shows data for a date specified in its title (YYYYMMDD).

Figure 6 .
Figure 6.The relationship between burn severity categories and dLST between pre-fire (13 July 2009) and post-fire (29 July 2009) images.Bars indicate the confidence interval of the average values (α = 0.01).

Figure 10 .
Figure 10.Mean LST values in different burn severity categories by date and aspect intervals (degrees).The title of each graphic indicates the image date (YYYYMMDD).

Atlantic Ocean Me dit err an ea n Se a Atlantic Ocean Me dit err an ea n Se a
RH) obtained from the Hurdes-Azabal meteorological station[57].The station is part of the Spanish Agroclimatic Information System for Irrigation (SIAR)[58]and is about 10 km from the study site.

Table 1 .
Landsat-5 TM images and meteorological data on the dates involved.RH, relative humidity.

Table 3 .
Average NDVI values by fire severity category and date (MMDD).