A Framework for Multi-Dimensional Assessment of Wildfire Disturbance Severity from Remotely Sensed Ecosystem Functioning Attributes

Abstract

Wildfire disturbances can cause modifications in different dimensions of ecosystem functioning, i.e., the flows of matter and energy. There is an increasing need for methods to assess such changes, as functional approaches offer advantages over those focused solely on structural or compositional attributes. In this regard, remote sensing can support indicators for estimating a wide variety of effects of fire on ecosystem functioning, beyond burn severity assessment. These indicators can be described using intra-annual metrics of quantity, seasonality, and timing, called Ecosystem Functioning Attributes (EFAs). Here, we propose a satellite-based framework to evaluate the impacts, at short to medium term (i.e., from the year of fire to the second year after), of wildfires on four dimensions of ecosystem functioning: (i) primary productivity, (ii) vegetation water content, (iii) albedo, and (iv) sensible heat. We illustrated our approach by comparing inter-annual anomalies in satellite-based EFAs in the northwest of the Iberian Peninsula, from 2000 to 2018. Random Forest models were used to assess the ability of EFAs to discriminate burned vs. unburned areas and to rank the predictive importance of EFAs. Together with effect sizes, this ranking was used to select a parsimonious set of indicators for analyzing the main effects of wildfire disturbances on ecosystem functioning, for both the whole study area (i.e., regional scale), as well as for four selected burned patches with different environmental conditions (i.e., local scale). With both high accuracies (area under the receiver operating characteristic curve (AUC) > 0.98) and effect sizes (Cohen’s |d| > 0.8), we found important effects on all four dimensions, especially on primary productivity and sensible heat, with the best performance for quantity metrics. Different spatiotemporal patterns of wildfire severity across the selected burned patches for different dimensions further highlighted the importance of considering the multi-dimensional effects of wildfire disturbances on key aspects of ecosystem functioning at different timeframes, which allowed us to diagnose both abrupt and lagged effects. Finally, we discuss the applicability as well as the potential advantages of the proposed approach for more comprehensive assessments of fire severity.

1. Introduction

Wildfire disturbances are considered an integral part of the natural dynamics of ecosystems in several biomes [1,2]. Notwithstanding, wildfire events can modify the composition, structure, and functioning of ecosystems and, consequently, the provision of ecosystem services to humankind [3]. Among those alterations, changes in ecosystem functioning are of particular interest, as fire can cause rapid modifications in the matter and energy budgets of ecosystems [4]. This is because attributes of ecosystem functioning have a shorter response time to disturbances than structural or compositional attributes of landscapes and are more directly connected to ecosystem services [5].

Wildfires can have profound impacts on many key aspects of the flows of matter and energy. These include dimensions of ecosystem functioning related to the biogeochemical cycles of carbon (e.g., primary productivity, biomass), and water (e.g., vegetation water content, soil moisture), as well as energy balances (e.g., albedo, latent heat, sensible heat). Namely, fires play an important role in the terrestrial biosphere carbon cycle [6], as they can cause substantial losses in carbon storage [7], including both aboveground biomass [8] as well as soil organic matter [9]. Moreover, fire-mediated changes in nutrient concentrations can ultimately limit productivity [10] as well as induce phenological changes [11]. Wildfires can also affect both water supply [12] and quality [13,14], as well as induce vegetation mortality due to prolonged loss of foliar moisture [15]. Additionally, wildfires can have negative impacts on evapotranspiration [16,17], which, in turn, also influences radiation exchanges, as evapotranspiration is related to latent heat flux [18]. Changes in wildfire frequency and/or severity may also have strong impacts on Earth’s surface radiative budget [19,20]. Burn severity influences the magnitude of changes in land surface albedo [21,22], which can undergo a darkening and/or a brightening that can persist in time [23,24]. Such alterations cause an increase in land surface temperature immediately following fires [25], which influences both burned area and the duration of fires [26], since vegetation is a regulator of land surface energy fluxes [27]. Although the effects of wildfires span across multiple dimensions of ecosystem functioning, most studies evaluating those effects seldom address those multiple dimensions at the same time, thus failing to fully depict spatiotemporal changes caused by these disturbances and limiting the assessment of environmental impacts.

As wildfire events have been increasing in recent decades, both in terms of frequency and intensity [28,29], exacerbated by global climate change and by shifts in land use and forest management [30], there is an increasing need for methods to assess and monitor the ecological consequences of such disturbances [31]. Over the last decades, remote sensing (RS) has revolutionized the way environmental changes are monitored, with new satellite sensors (e.g., Sentinel-2 MSI, Landsat-8 OLI/TIRS) providing valuable data at increasingly higher temporal, spatial, and spectral resolutions. Other satellites provide long-term archives of images that are useful to support temporal studies (e.g., Terra/Aqua MODIS, Landsat 5/7 missions). This growing variety and quality of remotely sensed data enables a wide range of fire-related applications, from the detection of fire occurrences [32] to an evaluation of fire severity [33] and monitoring of post-fire recovery [31,34,35] and resilience [36]. The ability to assess and map the heterogeneous (both spatially and temporally) effects of wildfires on ecosystems with remotely sensed data is a major asset not only for risk assessment and governance but also for post-fire management and restoration [30,37,38,39]. In this regard, spectral indices and metrics derived from satellite multi-spectral images have been used to assess different aspects of fire severity, usually through the comparison of pre-fire versus post-fire values. For instance, spectral indices sensitive to photosynthetically active radiation and/or water content, such as the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR) have been used, in combination with data collected in-field (e.g., through the Composite Burn Index; CBI), to assess “burn severity” (e.g., [39,40,41]), which is a post-fire measure of organic matter loss and soil alteration [30,37,42]. Furthermore, both satellite-derived land surface temperature (LST; [43]) and albedo [22] have been used for the same end. Nevertheless, to fully characterize and understand the effects of wildfire disturbances on multiple dimensions of ecosystems, a better translation of spectral indices into informative ecosystem variables is needed. More diverse sets of indicators are thus required to extend “burn severity” assessments to a wider variety of fire effects, i.e., towards more comprehensive approaches in “fire severity” assessment.

In recent decades, Ecosystem Functioning Attributes (EFAs) derived from satellite image time-series (SITS) have been increasingly used in a wide range of ecological applications due to their strong relation to biophysical properties and processes of ecosystems [5]. These include wildfire-related applications, such as improving the detection of wildfire disturbances in space and time [44], but also a wide range of other ecological applications, such as describing major ecological patterns [45,46,47], predicting and projecting species distributions [48,49,50], and supporting the definition of conservation priorities [51]. More specifically, EFAs extracted from SITS can provide information on inter-annual quantity (e.g., mean, minimum, maximum), seasonality (i.e., seasonal range or variability), and timing (e.g., dates of specific moments) components over multiple dimensions of ecosystem functioning. These include frequently used satellite-derived descriptors of the intra-annual dynamics of carbon gains, such as primary productivity, vegetation seasonality, and phenology (e.g., [46,52,53]). However, analogous measures can also be extracted from RS-derived proxy descriptors related to other dimensions of ecosystem functioning, such as water content, albedo, and sensible heat (e.g., [50,54]). By combining these four complementary dimensions of ecosystem functioning across each component (i.e., quantity, seasonality, and timing), remotely sensed EFAs can enable in-depth and integrative characterizations of fire severity [3,44].

The main objective of this study is to propose, describe, and showcase a framework to support enhanced assessments of fire severity, encompassing multiple dimensions of ecosystem functioning. Our approach is based on the effects of wildfire disturbances, at short (i.e., the year of the fire) to medium term (i.e., up to the second year after the fire), on descriptors of the intra-annual dynamics (i.e., EFAs) of four essential dimensions of the flows of matter and energy in ecosystems: (i) primary productivity; (ii) vegetation water content; (iii) albedo; and (iv) sensible heat. To illustrate our approach, we assessed the potential of inter-annual anomalies (i.e., deviations from the normal inter-annual variability) of a large set of EFAs extracted from MODIS data, for the period between 2000 and 2018 in the NW Iberian Peninsula, as estimators of the effects of wildfire disturbances. To this end, we compared, ranked, and analyzed the main patterns of those inter-annual anomalies at the regional scale. We then showcased the translation of inter-annual anomalies of the selected EFAs into indicators of wildfire disturbance severity for four selected burned patches within the study area and compared and analyzed their main spatial and temporal patterns at the local scale. Finally, we discussed the potential and added value of the proposed approach to improve RS-based assessment, mapping, and monitoring of wildfire disturbance severity on multiple dimensions of ecosystem functioning.

2. Materials and Methods

2.1. Generic Framework

2.1.1. General Workflow

The proposed approach to evaluate the effects of wildfire disturbances on ecosystem functioning consists of using indicators derived from Ecosystem Functioning Attributes (EFAs) for each of the four key dimensions of ecosystem functioning (i.e., primary productivity, vegetation water content, albedo, and sensible heat). Figure 1 presents the general workflow of the proposed approach, which is composed of five steps, described in the following sub-sections.

2.1.2. Step 1—Satellite Time-Series

First, satellite image time-series (SITS) for each of the four dimensions of ecosystem functioning have to be collected from available satellite products and preprocessed (Figure 1—”Step 1”). From the large number of spectral band combinations and vegetation indices available that have been successfully used for fire applications (e.g., [43,44,55,56,57]), we propose, for that end, the use of land surface temperature (LST) for the “heat” dimension, as well as the Tasseled Cap Transformation (TCT) features of “Brightness” (TCTB), “Greenness” (TCTG), and “Wetness” (TCTW) for the “albedo”, “primary productivity”, and “vegetation water content” dimensions, respectively. LST is a measure of the thermal energy emitted by different surfaces, calibrated from thermal emissivity measured by the satellite instrument [58]. On the other hand, TCT features are well-known, sensor-specific, linear combinations of bands in the visible, near-infrared, and short-wave infrared regions of the electromagnetic spectrum [59]. These three TCT features of “Brightness”, “Greenness”, and “Wetness” result from a rotation of principal component axes, derived from a global sample, to be aligned with the biophysical parameters of albedo, amount of photosynthetically active vegetation, and soil moisture, respectively [55].

2.1.3. Step 2—Extraction of EFAs

Next, metrics describing aspects of the intra-annual dynamics of ecosystem functioning (i.e., EFAs) have to be extracted from the SITS for each target year (Figure 1—”Step 2”). These types of metrics are commonly used to describe intra-annual aspects of ecosystem functioning related to quantity, seasonality, and timing (e.g., [5,46,48]). Examples of “quantity” metrics are the annual mean or median, the maximum or minimum values, or an integrated cumulative sum during a particular part of the year (e.g., the growing season). For “seasonality”, metrics such as the intra-annual standard deviation or range can be used. Examples of “timing” metrics are the dates of the annual maximum, minimum, or other important moments (e.g., the start/end of the growing season).

2.1.4. Step 3—Computation of EFA Anomalies

To obtain inter-annual EFA anomalies, it is necessary to first compute reference values for each pixel (Figure 1—”Step 3a”). To this end, EFAs are extracted for a multi-annual reference profile, which is obtained by averaging all years, after excluding values for years with any record of wildfire occurrences (this procedure is similar to the one used in [34,60], except for the exclusion of values in years with fire occurrences). Then, the anomalies are calculated by subtracting the reference values to those of each year, for each EFA (Figure 1—”Step 3b”). More precisely, in the case of linear (i.e., non-circular) metrics, the anomalies are calculated using the following expression:

$$A_m=X_m-R_m,$$

(1)

where A_m represents the annual anomalies of EFA m; X_m is the annual values of EFA m, and R_m represents the values of the reference year for EFA m. In the case of circular metrics—e.g., a day of the year (DOY)—anomalies are calculated through the smallest differences between angles, in the following manner:

$$A_m=\left(\frac{365}{2\pi }\times \theta \right)+35,$$

(2)

where A_m represents the annual anomalies of EFA m, and theta is defined as follows:

$$\theta =\{\begin{array}{cc}arctan\left(\frac{y}{x}\right),& ifx>0\\ arctan\left(\frac{y}{x}\right)+\pi ,& ifx<0\mathrm{and}y⩾0\\ arctan\left(\frac{y}{x}\right)-\pi ,& ifx<0\mathrm{and}y<0\\ \frac{\pi }{2},& ifx=0\mathrm{and}y>0\\ \frac{-\pi }{2},& ifx=0\mathrm{and}y<0\\ undefined,& ifx=0\mathrm{and}y=0\end{array}$$

(3)

with:

$$x=cos\left(X_m-R_m\right),$$

(4)

and:

$$y=sin\left(X_m-R_m\right),$$

(5)

where X_m and R_m represent the same as in Equation (1).

2.1.5. Step 4—EFA Ranking and Selection Procedures

Since multiple aspects of the intra-annual dynamics of ecosystem functioning can be obtained through the computation of EFAs, a selection process is necessary to reduce the initial set of metrics to an essential (i.e., relevant and less correlated) parsimonious set. In this selection process, classification models are used to rank all the EFAs considered in terms of their potential to measure the remotely sensed effects (as proxies of severity) of wildfire disturbances on the different EFAs (Figure 1—”Step 4a”). To this end, burned-area maps considered the response variable, while the anomalies of all EFAs are considered the predictive variables. The resulting ranking of EFA anomalies is used to compare all EFAs across (i) base SITS, with each corresponding to a different dimension of ecosystem functioning (i.e., “primary productivity”, “vegetation water content”, “albedo”, and “sensible heat”); and (ii) the specific metric (e.g., mean vs. median) and (iii) component of the inter-annual dynamics of ecosystem functioning (i.e., “quantity”, “seasonality”, and “timing”) used, each corresponding to an EFA. Furthermore, potential temporally lagged effects can also be evaluated by including values of EFA anomalies of different years concerning the year of the wildfire disturbance event (e.g., the anomalies for the year of the fire event, plus the first year after, the second year after, etc.).

The resulting ranking can be coupled with additional selection criteria, such as pairwise correlations, to minimize potential redundancy between EFAs so that a compact set of twelve EFAs (i.e., one for each dimension and each component of the inter-annual dynamics of ecosystem functioning) are selected (Figure 1—”Step 4b”). Then, both the magnitude and the direction of the pairwise relationship between the response variable (e.g., burned areas) and each of the selected EFAs are assessed (e.g., using Cohen’s d effect size statistics [61]). This allows for the inter-comparison of selected EFAs across all the dimensions and components of ecosystem functioning, as well as other factors considered (e.g., the year of the early post-fire period). Based on these procedures, a final selection of EFAs is carried out to support the analysis of the main patterns of wildfire disturbance severity on the dimensions of ecosystem functioning considered.

2.1.6. Step 5—Translation into Indicators of Wildfire Disturbance Severity

To support a multi-dimensional assessment of the effects of wildfire disturbances on ecosystem functioning, the final set of selected EFAs are then converted into indicators of wildfire disturbance severity (Figure 1—”Step 5”). For obtaining these indicators, the absolute values of the corresponding EFA anomalies are used, as these constitute measures of the magnitude of deviations from reference, thus allowing easier comparison across different dimensions. The resulting compact set of indicators is then used to support more detailed analyses.

2.2. Test Case

2.2.1. Study Area

To illustrate our proposed approach, we chose the northwest of the Iberian Peninsula (NW-IP) as the overall study area (Figure 2), since this is a hotspot in terms of wildfire occurrences within the European context, being a highly fire-prone landscape [3]. This area extends to the north and the west to the Atlantic Ocean, while its southern boundary was delimited by the edge of MODIS tile h17v04. Its eastern boundary, on the other hand, was defined to include the northern half of mainland Portugal, as well as the Spanish autonomous community of Galicia.

The NW-IP features strong environmental gradients and includes a major biogeographic transition, from an Atlantic climate with temperate deciduous broadleaf and mixed forests in the north and northwest, to a Mediterranean climate with evergreen sclerophyllous vegetation towards the southeast. Land cover and uses are well diversified and highly heterogeneous (

Table 1. Summary of baseline environmental characteristics of the overall study area for regional-scale analyses (northwest Iberian Peninsula, NW-IP), as well as the four selected burned patches selected for local-scale analyses (A–D). Elevation was extracted from MERIT DEM [64]; climate variables were extracted from the CHELSA Bioclim dataset [65]; the percentages of land-cover classes were extracted from CORINE Land Cover (CLC) 2000 [66].

	NW-IP	A	B	C	D
Year of fire	―	2003	2005	2005	2013
Average burned area (ha)	3617	14,625	17,600	19,325	14,850
Distance to coast (km)	―	212	9	91	153
Average elevation (m a.s.l.)	580	801	261	712	505
Average temperature (°C)	12.7	13.0	14.1	12.5	14.1
Minimum temperature (°C)	3.1	1.2	7.5	1.0	2.8
Maximum temperature (°C)	2.5	28.9	21.8	26.7	29.1
Average total precipitation (mm·yr−1)	1139	1075	1747	1229	620
% of Urban	1.8	0.0	1.1	0.1	0.2
% of Agricultural	29.0	9.3	8.0	5.8	26.0
% of Broad-leaf forests	7.1	0.8	1.6	0.2	4.9
% of Coniferous forests	5.9	34.5	17.7	32.8	10.8
% of Mixed forests	8.7	4.5	9.1	3.4	1.5
% of Natural grasslands	5.6	0.2	4.0	0.0	17.4
% of Shrublands	20.2	50.1	32.7	57.1	38.2
% of Bare rocks or sparsely vegetated	2.0	0.6	25.7	0.6	1.0

), including cropland areas, urban areas, plantation forests (mostly maritime pine, eucalyptus, and mixed stands), and shrublands in different successional stages, along with historical use of fire for agrosilvopastoral purposes. Wildfire occurrences are more frequent during the months between July and October—but especially in August—with burned areas over 10,000 ha being increasingly common. Wildfire events have been increasing in recent decades because of changes in land use, which mainly resulted from extensive abandonment of farming and husbandry, increasing the fuel load and its continuity in the landscape [30,62]. The frequency of extreme fire events, also known as mega-fires (i.e., complex catastrophic events in terms of size, intensity, resistance to control, severity, etc. [30]), has also been increasing, despite the huge investments in fire suppression [28], making this one of the regions with both the highest annual values of burnt area in Europe and the highest density of ignitions among southern European countries [63].

To showcase our proposed framework, the full study area in NW-IP, including both burned and unburned areas of NW-IP between 2000 and 2018, was used for regional-scale analyses. Then, for local-scale analyses, four individual burned patches (see boxes labelled A–D in Figure 2) were selected from different geographical locations within the study area, featuring diverse baseline environmental characteristics (e.g., altitude, distance to coast, climate, proportions of major vegetation types), and with the focal wildfire event having occurred at different years within the considered period (A: 2003; B and C: 2005; D: 2013). Table 1 presents a summary of the characteristics of both the full study area, used for regional-scale analyses (i.e., NW-IP), as well as the four individual burned patches used for local-scale analyses.

2.2.2. Satellite Data Preprocessing

In our test case, we used SITS extracted from MODIS products. Most tasks were undertaken within the R statistical programming environment [67], using mainly the “raster” package [68], with additional R packages complemented by other software for more particular tasks, as specified whenever relevant.

We started by downloading and preprocessing all available images encompassing NW-IP (i.e., MODIS tile h17v04), between 2000 and 2018, from the following MODIS products, using the “MODIStsp” R package [69]:

• MODIS/Terra+Aqua Burned Area Monthly L3 Global 500m (MCD64A1), Collection 6 [70];
• MODIS/Terra Land Surface Temperature/Emissivity 8-Day L3 Global 1km (MOD11A2), Collection 6 [71];
• MODIS/Terra Surface Reflectance 8-Day L3 Global 500m (MOD09A1), Collection 6 [72].

Preprocessing consisted of converting to the GeoTIFF file format, re-projecting to a common reference system (WGS84/UTM29N), and resampling to 500-m pixel size using the nearest neighbor interpolation algorithm. Additionally, a filter based on the Hampel outlier identifier [73,74] with a window size of 5 was applied to the resulting SITS to reduce residual noise that hinders time-series data.

Next, LST and TCTB, TCTG, and TCTW were extracted from the MOD11A2 and MOD01A1 products, respectively, using “GDAL” [75] and the “rasterio” Python package [76]. For LST, only the day temperatures were used, which were then rescaled to degrees Celsius. The three TCT features of “Brightness”, “Greenness”, and “Wetness” were obtained by combining the 7 bands available in the MOD09A1 product, using the coefficients presented in

Table 2. MODIS-specific coefficients to calculate the Tasseled Cap Transformation (TCT) features “Brightness”, “Greenness”, and “Wetness”. These features result from a rigid rotation of principal component axes, which are aligned with the biophysical parameters of albedo, the amount of photosynthetically active vegetation, and soil moisture, respectively [59].

Band			Coefficients
No.	Name	Range (nm)	Brightness	Greenness	Wetness
1	Red	620–670	0.4395	–0.4064	0.1147
2	NIR 1	841–876	0.5945	0.5129	0.2489
3	Blue	459–479	0.2460	–0.2744	0.2408
4	Green	545–565	0.3918	–0.2893	0.3132
5	NIR 2	1230–1250	0.3506	0.4882	–0.3122
6	SWIR 1	1628–1652	0.2136	–0.0036	–0.6416
7	SWIR 2	2105–2155	0.2678	–0.4169	–0.5087

NIR: near-infrared; SWIR: short-wavelength infrared.

[59].

Burned areas extracted from the MCD64A1 product were filtered by size, eliminating all burned patches smaller than 100 ha, following a set of criteria detailed in previous work (see [44]). To further reduce the remaining noise in the data, a Whittaker–Henderson smoother [77,78] was applied to the four extracted SITS, with a lambda parameter equal to 2 (a rather low value, to conserve the characteristics of the time-series while still reducing extreme noise).

2.2.3. EFA Anomalies Computation

A total of 60 EFAs were extracted, using 15 different intra-annual metrics (

Table 3. List of metrics extracted from each of the four satellite image time-series used in this study (i.e., land surface temperature and the Tasseled Cap Transformation features of “Brightness”, “Greenness”, and “Wetness”); as well as the number of predictive variables for models extracted from each one.

Component	Metric	Abbreviation	No. of Variables for Models
Quantity	Mean (or average)	avg	12
	Median	mdn	12
	Maximum	max	12
	Minimum	min	12
Seasonality	Standard deviation	std	12
	Median absolute deviation	mad	12
	Absolute range	rng	12
	Relative range	rrl	12
	Non-parametric relative range	rnp	12
Timing	Time (of the year) of maximum	tmx	12
	“Winterness” of maximum	wmx	12
	“Springness” of maximum	smx	12
	Time (of the year) of minimum	tmn	12
	“Winterness” of minimum	wmn	12
	“Springness” of minimum	smn	12
Total			180

), from each of the four SITS—LST, TCTB, TCTG, and TCTW—for each year between 2000 and 2018 for NW-IP. Among these, well-known measures were computed for “quantity”—mean, median, maximum, and minimum annual values—and “seasonality”—standard deviation, median absolute deviation, and absolute range. In addition to those, the relative range and a non-parametric relative range were also computed using the following expressions:

$$rrl=lo{g}_{10}\left(\left|\frac{rng}{avg}\right|\right)$$

(6)

and

$$npr=lo{g}_{10}\left(\left|\frac{rng}{mdn}\right|\right),$$

(7)

where rrl is the relative range metric, while rng is the absolute range, avg is the mean, npr is the non-parametric relative range, and mdn is the median.

As for “timing”, the times (of the year) of the maximum and minimum values were used. However, as variables of time of the year are circular by definition, we also tested two linearized versions for each, called “winterness” and “springness”, using the following expressions:

$$\mathrm{winterness}=cos\left(\frac{2\pi }{365}\times \left(x-35\right)\right)$$

(8)

and:

$$\mathrm{springness}=sin\left(\frac{2\pi }{365}\times \left(x-35\right)\right),$$

(9)

where x is the day of the year (DOY) of either the maximum or minimum value (tmx or tmn, respectively). Note that these last two expressions (i.e., Equations (8) and (9)) include the conversion from DOY to radians, assuming that a full year has 365 days; as well as a rotational realignment of -35 days so that the maximum and minimum values of winterness and springness (i.e., 1 and −1, respectively) correspond to the mid-point of each season (in the case of NW-IP, as it is located in the northern hemisphere, the mid-point of winter is in early February). These metrics can be interpreted as corresponding to the temporal proximity, within the year, to either the start of winter or spring, respectively.

Finally, anomalies were then obtained for all EFAs for each of the four SITS by calculating either linear differences (Equation (1)), for all metrics except “tmx” and “tmn”, or the smallest differences between angles, for those two circular metrics (Equations (2)–(5)).

2.2.4. Ranking and Selection of EFAs

A ranking of EFA anomalies was obtained for NW-IP through a classification modeling procedure, in which the response variable was the binary burned class (i.e., burned vs. unburned) extracted from the MODIS burned area product. As predictive variables for these models, we used the inter-annual anomalies in each one of the 60 EFAs (i.e., each of the 15 intra-annual metrics described in Table 3, for each of the four SITS), for each one of three years since the fire occurrence: the year of the fire, as well as the first and second years after (i.e., years 0, +1, and +2). This was to allow for the evaluation of potential temporally lagged effects (for which only the years from 2003 to 2016 were used). In total, the number of predictive variables used in these models was 180 (see Supplementary Information S1 for additional information on the predictive variables used in these models). To this end, the Random Forest (RF) “ensemble learning” algorithm was used, since this technique scales well on large and/or high-dimensional data [79], thus allowing for the full set of variables to be tested without the need for a selection before modeling.

As the number of pixels identified as “burned” represented only ca. 1.3% of the total overall study area (i.e., NW-IP), data balancing through down-sampling of the majority class (i.e., “unburned”) was applied. To that end, 100 random sample sets were defined, consisting of all the pixels identified as “burned”, as well as an equal number of randomly selected “unburned” pixels for each year, in a total of 122,658 observations (i.e., n = 61,329, per class). Then, an RF classification model was calibrated for each one of the 100 sample sets, using the “caret” R package [80], with the hyperparameters set to the default values (including a fixed value for “mtry” equal to the default squared root of the number of predictive variables), and under a leave-group-out cross-validation (LGOCV) scheme with 50:50 train/test split in each one of 10 repetitions.

Model performance was then evaluated using measures specific to two-class problems—sensitivity (or true positive rate), specificity (or true negative rate), and the area under the receiver operating characteristic curve (AUC)—which were then averaged across the 100 resulting models. Finally, variable importance scores—a common approach to extract interpretable information on the contribution of different variables from RF models [79]—based on the mean decrease in accuracy, were extracted and also aggregated to the median values across the 100 sample sets to support the EFA ranking.

Based on the ranking of variable importance scores obtained from the RF models, one EFA was selected for each dimension and each component, giving a total of 12, to support the analysis of the main patterns of wildfire disturbance severity on multiple dimensions of ecosystem functioning in NW-IP.

2.2.5. Analysis of Indicators of Wildfire Disturbance Severity

To support the analysis of both the directionality (i.e., increasing or decreasing) and magnitude of the pairwise relationships between the regional-scale burned areas and each of the selected EFAs, Cohen’s d effect size statistics [61] and Spearman pairwise rank correlations were calculated for the 100 sample sets. From this, as well as from the distributions of the selected EFA anomalies, the main overall patterns of the effects of wildfire disturbances on ecosystem functioning were analyzed for NW-IP and compared across the four dimensions and three components of ecosystem functioning, as well as across years of the early post-fire period.

It should be noted that effect sizes were applied here to avoid the use of statistical hypothesis testing based on significance (i.e., p-values), which have important known limitations such as not measuring the effect size or the importance of a result [81] or the strong influence of sample size, which is hugely influential in determining significance levels [82] (i.e., when sample sizes are large enough, any effect can produce a small p-value if the sample size or measurement precision is high enough [81]; with large sample sizes, often happening at a pixel-level analysis of SITS).

We then computed four indicators of the effects of wildfire disturbances on ecosystem functioning (i.e., one for each dimension) for four selected burned areas (see Figure 2) at a local scale to illustrate the proposed approach for RS-based assessment, mapping, and monitoring of wildfire disturbance severity on multiple dimensions of ecosystem functioning. The selection of EFAs for these final indicators took into consideration the results from both the Random Forest model-based ranking and the pairwise correlations between EFAs. The translation of the selected EFAs into indicators was done by taking the absolute values of the anomalies of the selected EFAs for the four burned areas. Finally, spatial patterns between all indicators, across years, were further analyzed, local-scale, through local Spearman correlations, with a neighborhood window size of 5 pixels.

3. Results

3.1. EFA Ranking

In terms of the performance of the Random Forest (RF) models for our test case, the results obtained for leave-group-out cross-validation can be considered very good across all sample sets, with around 0.9815 ± 0.0004 for AUC, 0.9080 ± 0.0015 for sensitivity, and 0.9647 ± 0.0010 for specificity.

As for variable importance scores, EFAs extracted from TCTG had overall higher values than the other dimensions (Figure 3; see Supplementary Information S2 for the full plot of variable importance scores for all non-aggregated variables). Considering each of the three different EFA components, the ones measuring aspects of “quantity”—particularly the mean and either one of the extreme values (minimum or the maximum, depending on the specific case)—obtained overall higher variable importance scores than “seasonality” or “timing” EFAs, with the notable exception of seasonality EFAs extracted from TCTG at year 0. Conversely, the same exception was observed within all seasonality EFAs, where otherwise, LST was the dimension from which EFA anomalies obtained higher variable importance scores overall, especially for standard deviation. As for timing EFAs, the scores obtained were generally lower than for the two other groups, except for the ones extracted from either TCTG or LST at year 0, particularly the springness of the time of the year of the minimum TCTG value and the winterness of the time of the year of the maximum LST value, closely following the corresponding quantity EFAs.

Across the four different dimensions of ecosystem functioning, those extracted for the year of fire (i.e., year 0) generally obtained better scores than those of the two subsequent years. However, for quantity EFAs extracted from TCTW, the first year after the fire (i.e., year +1) obtained better scores. Furthermore, quantity EFAs extracted from TCTB also achieved higher scores for the second year after the fire (i.e., year +2), although not as high as for year 0. Finally, in terms of individual EFAs, the top-ranked ones were extracted from TCTG, with the minimum value (TCTG-min) and the standard deviation of TCTG (TCTG-std) at year 0 holding the highest scores among all EFAs tested (see Supplementary Information S1 for additional information on the ranking results).

3.2. Analysis of Effects

Considering the RF-based importance ranking, 12 EFAs (

Table 4. List of selected EFAs with the corresponding effects categories, based on Cohen’s d effect size statistic. The number of arrows in the “Effect category” columns symbolizes the magnitude of the effect size, with |d| < 0.2 “negligible”, |d| < 0.5 “small”, |d| < 0.8 “medium”, and otherwise “large”. The direction and color of the arrows symbolize the sign of the effect size, with ascending arrows depicting positive (non-“negligible”) effect sizes while descending arrows depict negative (non-“negligible”) effect sizes (“negligible” effects are symbolized by a horizontal dash).

Dimension	Component	Attributes (EFA)	Effect Category
			Year 0	Year +1	Year +2
Primary productivity	Quantity	TCTG-min	↘↘↘	↘↘↘	―
	Seasonality	TCTG-std	↗↗↗	↗↗	―
	Timing	TCTG-tmn	↘↘↘	―	―
Vegetation water content	Quantity	TCTW-avg		↘↘↘	↘↘↘
	Seasonality	TCTW-std	―	↗↗	―
	Timing	TCTW-tmn	―	―	―
Albedo	Quantity	TCTB-avg	↘↘↘	↗	↗↗↗
	Seasonality	TCTB-std	―	↗	―
	Timing	TCTB-tmx	―	―	―
Sensible heat	Quantity	LST-max	↗↗↗	↗↗↗	↗
	Seasonality	LST-std	↗↗	↗↗↗	↗
	Timing	LST-tmx	↗↗↗	―	―

) were selected, one for each dimension (i.e., “primary productivity”, “vegetation water content”, “albedo”, and “sensible heat”) and each component of intra-annual dynamics (i.e., “quantity”, “seasonality”, and “timing”). As for “timing” EFAs, it should be noted that the selected metrics do not correspond to the highest-ranked ones for each case. In this instance, parsimony and interpretability criteria were also considered, since linearized versions of timing EFAs (i.e., “winterness” and “springness”) often require the use of both at the same time to better understand in which month/season the minima/maxima occur (see Supplementary Information S3 for partial dependence plots of the predictive variables of the models, extracted from the 12 selected EFAs).

Table 4 summarizes the results for the effect sizes (both magnitudes and directions) between the response variable and each of the selected EFA anomalies (see Supplementary Information S4 for additional information on the effect size results). These showed overall strong negative effects (i.e., lower values associated with wildfire disturbances) for TCTG-min and TCTW-avg (especially at years 0 and +1 and years +1 and +2, respectively), while positive (i.e., raised values associated with wildfire disturbances) for LST-max (especially at years 0 and +1). For TCTB-avg, however, the effects changed from strongly negative at year 0 to strongly positive at year +2. As for seasonality EFAs, the effects were, overall, positive, with the strongest magnitude associated with TCTG-std and LST-std. Finally, timing EFAs had overall weaker effects, except for TCTG-tmn (with large negative effects) and LST-tmx (with large positive effects), both only at year 0.

3.3. Main Patterns in EFA Anomalies

The distributions of the deviations from the normal variation (i.e., anomalies) for each of the twelve selected EFAs, grouped by year after fire (i.e., years 0, +1, and +2) and by burned vs. unburned pixels, are shown in Figure 4 (see also Supplementary Information S5 for the distributions of the reference values of all pixels in the NW-IP).

Overall, there were notable differences in the distributions obtained for burned vs. unburned pixels, with the most notable contrast having been observed for TCTG-tmn at year 0, with median anomalies corresponding to occurrences of 80–88 days earlier than the reference. Within the same dimension, differences were also notable for TCTG-min (especially at years 0 and +1), and also for TCTG-std (especially at year 0), with values of burned pixels considerably below/above (respectively) those of unburned pixels. Considerable differences between burned vs. unburned were also obtained for TCTW-avg, but only at years +1 and +2, while for TCTB-avg, they were observable at years 0 and +2, with opposite directions (i.e., negative vs. positive). As for heat-related EFAs, differences in the distributions of anomalies between burned and unburned pixels were notable for LST-max and LST-std across all years after fire, but especially at years 0 and +1. Finally, LST-tmx showed a contrast between burned and unburned pixels, with median anomalies corresponding to delays of 8–16 days.

3.4. Multi-Dimensional Assessment of Wildfire Disturbance Severity

For the RS-based assessment of wildfire disturbance severity on the four dimensions of ecosystem functioning considered, at the local level, one indicator for each dimension was used, based on the previously selected anomalies in quantity EFAs (i.e., TCTG-min, TCTW-avg, TCTB-avg, and LST-max). This selection took into consideration the results from both the model-based ranking (i.e., quantity EFAs ranked better overall than those of seasonality or timing) and pairwise correlations (up to ±0.83; see Supplementary Information S6).

Figure 5 shows the median (±0.5 median absolute deviations) profiles of each of the four selected indicators, across the three years after fire, for each of the four burned areas of interest (see Supplementary Information S7 for additional information on the spatial and temporal patterns of the indicators of wildfire disturbance severity in the four burned areas A–D). This figure illustrates different post-fire trajectories, with varying relations between the four dimensions of ecosystem functioning across individual burned areas. Specifically, different areas had different ratios between “Productivity” and “Heat” (e.g., A vs. D). Furthermore, “Water” seemed to be less affected in area D than in the other areas, while “Albedo” in area D exhibited a different post-fire trajectory relative to the other areas, with the value for year +1 being slightly higher than for year 0.

Finally, Figure 6 shows the variation in the spatial and temporal patterns of each of the four selected indicators of primary productivity, vegetation water content, albedo, and sensible heat, for each of the four individual burned areas of interest (A–D). Overall, some wildfire disturbance severity “hotspots” are observable, while not necessarily coinciding across dimensions (see also Supplementary Information S7). Moreover, different burned areas exhibited different patterns across time between dimensions of ecosystem functioning, with temporal trends being divergent in some cases while convergent in others (e.g., “Water” and “Albedo” dimensions in the areas of interest C vs. D), or with different relative preponderances (e.g., “Productivity” vs. “Heat” in the areas of interest A vs. D), further illustrating the patterns observed in Figure 5.

4. Discussion

4.1. Fire Severity Patterns in the NW Iberian Peninsula

4.1.1. Effects of Wildfires across Dimensions and Components

In our test case, we found that the associations between burned areas and inter-annual EFA anomalies were particularly strong for primary productivity and sensible heat, although albedo and vegetation water content also revealed important effects caused by wildfire disturbances. This is in line with previously published research, as the sudden removal of green vegetation due to fire usually translates into abrupt breaks that are commonly observable in time-series of spectral vegetation indices (e.g., TCTG). This removal can result in substantial carbon losses, both in terms of aboveground biomass [8] as well as soil organic matter [9,83]. Moreover, fire-mediated changes in nutrient concentrations can ultimately limit productivity over long periods [10]. Because vegetation is a regulator of land surface energy fluxes [27], the removal of vegetation causes an increase in observed LST induced by a reduction in evapotranspiration [25], thus increasing the sensible-to-latent heat ratio [17]. Furthermore, the observed effects in albedo are also consistent with studies reporting either a darkening or a brightening (or both) effect after a fire (e.g., [22,84]) due to the vegetation removal as well as the presence (or not) of black carbon in soot, which absorbs visible solar radiation [84]. Finally, effects on vegetation water content could be related to loss of moisture in canopy foliage due to damage caused by fire, eventually leading to vegetation mortality [15]. Indeed, our results show that using information from multiple dimensions of ecosystem functioning can yield improved results over the use of fewer—or only one—indices for enhanced detection, mapping, and evaluation of ecological effects of disturbances such as wildfires.

Across the three different components of ecosystem functioning considered, anomalies in quantity EFAs (e.g., mean, minimum, maximum) had overall stronger associations with wildfire disturbances than seasonality (e.g., standard deviation) or timing (e.g., time of minimum or maximum) ones. For severity assessment purposes, quantity and seasonality EFAs seemed to be somewhat redundant, as there can be, in some cases, high pairwise correlations between EFA anomalies of those two types. Conversely, for other types of disturbances such as land clearing for agriculture and ranging, the impact on seasonality (i.e., seasonal variability) can be much stronger than on quantities (see [85]). While some studies show that wildfire disturbances can have profound impacts on aspects of timing of key dimensions of ecosystem functioning (e.g., by inducing phenological changes [11]), our approach may not be best-suited to assess such modifications. This is because the intra-annual descriptors used (i.e., EFAs) seemed to more likely capture the location of extreme values (i.e., minimum/maximum) within the year corresponding to direct effects of the wildfire disturbance rather than translating changes (e.g., delays or lags) in the annual cycles of the underlying ecological processes.

As for the specific statistical measures used to extract EFAs, stronger effects can result from the use of parametric measures over that of their non-parametric equivalents. This may be because parametric measures are generally more sensitive (and susceptible) to extreme or abnormal values than non-parametric ones, which may or may not be desirable, depending on the specific application purpose. On the other hand, non-parametric measures may be more adequate for extracting long-term inter-annual trends than parametric ones (e.g., [5]).

The results obtained point to a strong added value of the proposed approach to discriminate the effects of wildfire disturbances on different aspects of ecosystem functioning.

4.1.2. Temporal Effects of Wildfires

Short-term effects of wildfire disturbances on ecosystem functioning seemed to attain better overall performance than medium-term effects, since inter-annual anomalies in intra-annual EFAs for the year of the fire occurrence obtained much better scores, overall, than those for the first and second years following fire. This is partially in line with other studies, as the majority of studies analyzing RS-based effects of fires focus on either short-term (i.e., abrupt) wildfire disturbance severity (e.g., [32,56]) or on longer-term effects more related to post-fire recovery (e.g., [3,34]), while the medium-term effects (i.e., from one or two years to five years after the fire) are often overlooked.

EFAs derived from TCTW (i.e., vegetation water content) seemed to be more affected by wildfire disturbances in the first year after the fire than the year in which the fire occurred, suggesting a temporal lag associated with this dimension relative to both the fire event and other dimensions. This observation (which, to our knowledge, has not been previously reported) could be linked to changes in hydrological dynamics, such as increased soil water repellency, precipitation run-off, and increased quantities of impervious materials such as ashes and soot [84], which can clog soil pores [86]. Alternatively, this effect could be due to post-fire changes in foliar moisture leading to mortality occurring over an extended period after fire [15], since TCTW is quite sensitive to the moisture content of vegetation [59].

As for the observed effects of wildfire disturbances on albedo (i.e., TCTB), this dimension of ecosystem functioning seemed to exhibit a shift in directionality, from negative in year 0 to positive in years +1 and +2, denoting a change in the relationships with the wildfire disturbance at short vs. medium term. This is also in line with the findings reported in other studies [22,87], which observed both a decrease in albedo (i.e., a darkening effect) immediately after the fire, and, usually, for a brief period, as well as a small increase (i.e., brightening) one year after the fire, which can be persistent. On the other hand, this darkening, which can be due to changes in the relative abundance of surfaces with distinct reflective properties [88] (e.g., ash, char, soot, bare soil), can, in some cases, persist for multiple years after a fire [23].

Finally, the effects of wildfire disturbances on both primary productivity and sensible heat (through TCTG and LST, respectively) seemed to be reflected across the three time periods analyzed. However, they were especially noticeable for the year of the fire occurrence. This also goes in the same direction as other studies addressing the impacts of fire on ecosystems, which report usually short-term effects for satellite-derived proxies of primary productivity (e.g., [41,89]) after fire. An increase in observed LST immediately following fires is also frequently reported in the literature [25]. This can sometimes still be observed one year following the fire [19], after which LST anomalies tend to become smaller, and seasonality starts governing the LST time-series [27].

Together, these results highlight the importance of taking into consideration both short- and medium-term effects of wildfire disturbances on multiple dimensions of ecosystem functioning.

4.1.3. General Patterns

Overall, EFA anomalies were able to efficiently discriminate burned and unburned areas, with very high values obtained for both effect sizes and performance measures of the Random Forest (RF) models between EFA anomalies (predictive variables) and burned areas (response variable). This allowed supporting the selection of EFAs from which a compact set of indicators were employed to assess the main patterns of wildfire disturbance severity in our test area.

Diverse spatio-temporal patterns of wildfire disturbance severity on the four dimensions of ecosystem functioning—primary productivity, vegetation water content, albedo, and sensible heat—could be observed both between as well as within individual burned areas across years (i.e., from the year of the fire to two years after). These results showcase the added value of (i) analyzing the effects of wildfire disturbances for multiple dimensions of ecosystem functioning, as opposed to addressing only overall burn severity in commonly used approaches (e.g., based on the Differenced Normalized Burn Ratio, ΔNBR); and (ii) taking into consideration both the short and medium term, thus enabling to obtain insights on potentially lagged effects of fires on ecosystem functioning, as different dimensions can exhibit different severities, and with different timings.

4.2. General Considerations about the Proposed Framework

4.2.1. Satellite Image Time-Series

While in the presented test case, we used data extracted from MODIS products exclusively, other sources of satellite image time-series (SITS) are available and can equally be applied using the proposed approach. Such data sources include platforms with diverse characteristics in terms of spatial, temporal, and spectral resolutions and length of the historical archive, such as SPOT-VEGETATION, PROBA-V, Landsat 5/7/8, or Sentinel-2/3 (e.g., [90]).

As for the base SITS to be extracted from those sources, we recommend, in the proposed approach, the use of LST and the TCT features of “Brightness”, “Greenness”, and “Wetness” as RS-based proxies of albedo, primary productivity, and vegetation water content, respectively. This option allows to derive information on the four above-mentioned dimensions of ecosystem functioning, from the smallest number of sources (the three TCT features are derived from the same source), compactly and coherently, and is often available for the major satellite data sources (e.g., [91]). However, in the case of primary productivity and vegetation water content, if data are not available to derive all four RS-based variables, alternative spectral indices could be used instead, such as the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI) for “productivity”, or the Normalized Difference Water Index (NDWI) or Land Surface Water Index (LSWI) for “water”.

4.2.2. Additional Data Sources

Well-designed field-based measurements could be used, if available, to validate the information provided by SITS, further enhancing the proposed approach, although it is not a mandatory requirement. To that end, data can be collected in the field, following robust sampling designs, such as spectral/radiometric readings, aerial surveys using cameras and sensors mounted on unmanned aerial vehicles (UAVs), or data collected under already existing and commonly used protocols (or adapted versions from these), such as the Composite Burn Index (CBI; e.g., [92]) or the improved geometrically structured Composite Burn Index (GeoCBI; e.g., [93]), to complement the information derived from SITS. Such data can, however, be difficult and/or costly to obtain, even more so if multiple observations are required to produce a coherent time-series. These types of field-level measurements are usually compared to commonly used indicators of burn severity (e.g., [40]) such as the differenced normalized burn ratio (ΔNBR)—or any of its enhanced forms, such as the relative normalized burn ratio (RdNBR) and the relativized burn ratio (RBR).

Unlike the more commonly used RS-derived methods to estimate burn severity, the approach proposed in this study can provide information on the effects of wildfire disturbances on multiple dimensions of ecosystem functioning—namely primary productivity, vegetation water content, albedo, and sensible heat—simultaneously, towards a comprehensive evaluation of fire severity.

4.2.3. Applicability and Future Directions

Multi-dimensional approaches to wildfire disturbances using EFA-like metrics are, to our knowledge, scarce (e.g., [94]). Very few studies have adopted a multi-dimensional approach using EFAs, regardless of the specific application (e.g., [48,50,95]). More studies should be conducted in the future, with diverse contexts in terms of baseline environmental conditions and fire disturbance regimes, to further test and improve the proposed approach. It is important to highlight that although the results obtained for the test case presented in this study may be somewhat specific to the study area and the analyzed period, the general approach could be applied, with relative ease, to other contexts and timeframes to conduct informed selections of the most efficient and informative EFA-based indicators for wildfire disturbance severity assessment.

The proposed approach makes use of annual descriptors of the dynamics of multiple dimensions of ecosystem functioning—so-called Ecosystem Functioning Attributes (EFAs)—extracted from SITS, to derive indicators of wildfire disturbance severity. However, another possible approach could be to derive indicators directly from the time-series, without the need to compute EFAs first. This would open the possibilities for other kinds of indicators, with finer temporal (i.e., <1 year) resolutions, on a (semi-)continuous basis, to be derived from SITS (e.g., [3]). On the other hand, although additional types of indicators could be obtained in this way, others would not—e.g., metrics of timing based on annual cycles, such as the date of the maximum or minimum, would only make sense on an annual basis. Besides, the annual nature of EFAs—and the derived indicators—facilitates both their computation as well as their interpretation, as it offers a better translation of spectral indices into informative ecosystem variables. This can be advantageous for some application purposes, such as for reporting from official authorities. Notwithstanding, approaches adopting either an annual or continuous basis could complement each other to characterize the effects of wildfire disturbances on ecosystem functioning.

Finally, in this study, we focused on short- to medium-term impacts of wildfire disturbances. However, long(er)-term effects of wildfire disturbance severity have been observed (e.g., [96]). Indeed, long-term impacts of wildfire disturbances on ecosystem functioning could also be monitored using the approach proposed in this study, provided that long-enough SITS are available. However, this was not tested in this study, as that would considerably increase both computational requirements as well as methodological complexity.

5. Conclusions

Here, we described a framework to assess the effects of wildfire disturbances on four key dimensions of ecosystem functioning—primary productivity, vegetation water content, albedo, and sensible heat. This approach is based on indicators derived from several descriptors of the intra-annual dynamics of ecosystem functioning—called Ecosystem Functioning Attributes (EFAs)—extracted from remotely sensed satellite image time-series (SITS). We found that all four dimensions of ecosystem functioning suffered important effects of wildfire disturbances—especially primary productivity and sensible heat. We also found that quantity EFAs held the highest potential for indicating wildfire severity. Finally, we found important differences in short- and medium-term effects between the four different dimensions of ecosystem functioning, suggesting temporally lagged effects in vegetation water content as well as directionality shifts in albedo. Together, these results highlight the added value of the proposed framework to enhance RS-based assessment, mapping, diagnostics, and monitoring of wildfire disturbances. Using this approach, a more comprehensive evaluation of fire severity can be attained, which is not captured by indicators derived from only one spectral index, measuring only overall burn severity. As such, this multi-dimensional framework can contribute to a deeper and more detailed understanding of the impact of fires in ecosystems, with the potential to support assessments of severity, recovery, and resilience and their interactions with biodiversity and ecosystem services.