**Preface**

Wildfires are one of the most significant disturbances worldwide, impacting both natural ecosystems and human communities. Many scientists and a portion of society believe that the frequency and intensity of fires are escalating due to shifts in land use, suppression policies, and climate change, leading to critical consequences for nature and the benefits it provides to people. Additionally, the expansion of human settlements into forested areas has altered our vulnerability to fire, further affecting the potential impacts of future events. In the context of these widespread changes, remote sensing emerges as a valuable tool with which to address the environmental and social challenges of forest fires and offer solutions. Its versatility, wealth of information, and rapid technological advancements make it indispensable for designing proactive strategies, monitoring risks, and analyzing damages over large areas, facilitating decision making and the development of pre- and post-fire management strategies.

This Special Issue aims to explore various aspects that contribute to the advancement of (I) fire driver characterization and the development of fire predictive models, (II) the assessment of burned areas (III) the assessment of burn severity and specific fire impacts, and (IV) the analysis of post-fire trajectories using remote sensing methods. The compilation comprises ten research articles addressing these four topics, utilizing a diverse range of methodologies and remote sensing platforms. The target audience for this Special Issue includes researchers, academicians, practitioners, and policymakers working in forest ecology, environmental science, disaster management, remote sensing, and geospatial technologies in general.

The success of this Special Issue can be attributed to the dedicated efforts and expertise of many individuals. First and foremost, we express our gratitude to the 45 authors for presenting their cutting-edge research and insights into remote sensing. Their exceptional contributions have enriched this Special Issue significantly. We would also like to extend our appreciation to the more than 100 reviewers who invested their time and knowledge in meticulously evaluating the submissions, thereby enhancing the quality of our Special Issue. Finally, we would like to acknowledge the editorial staff for their efficient management of the editorial processes, which facilitated the success of this Special Issue.

#### **V´ıctor Fern´andez-Garc´ıa, Leonor Calvo, Susana Suarez-Seoane, and Elena Marcos** *Editors*

## *Editorial* **Remote Sensing Advances in Fire Science: From Fire Predictors to Post-Fire Monitoring**

**Víctor Fernández-García 1,2,\*, Leonor Calvo 1, Susana Suárez-Seoane <sup>3</sup> and Elena Marcos <sup>1</sup>**


Fire activity has significant implications for ecological communities, biogeochemical cycles, climate, and human lives and assets. Approximately over half of the Earth's land surface is susceptible to fire, with around 3% experiencing annual burning according to coarse-resolution satellites [1], a value that is probably much higher according to recent estimates from finer satellite imagery [2]. Because of the vast extent of land burned over the world, landscape fires release approximately 23% of the global CO2 emitted annually from fossil fuels, modify Earth's energy fluxes through changes in surface albedo, and have an enormous influence on human health and the economy [1]. Fires also have a large influence on local ecosystems, affecting the ecosystem services provided to local communities. Thus, fires are a relevant phenomenon with an enormous area of impact every year.

Because of the relevance and impact of fires across the globe, the study of this phenomenon and the assessment of its consequences cannot be addressed only by field or laboratory studies. In this sense, the exploitation of remote sensing platforms, sensors, and methods is crucial to obtain accurate and extensive spatial and spatiotemporal information on fire and its consequences. Considering fire as a natural hazard, geo-spatial studies can be organized in multiple ways, one of them being based on the temporal point on which they focus in relation to fire. Thus, we can differentiate those studies focused on a pre-fire stage: for instance, those addressing topics that can help predict fire-related risks and thus are useful for pro-active management strategies. The second stage is the moment when fires occur. At that point, remote sensing might be useful to detect active fires, or to detect burn scars as evidence of fire. The next stage is the analysis of the immediate impacts and consequences of fire. This is related, but not limited to burn severity assessments, which indicate the overall environmental change caused by fire. The assessments immediately after fire are essential for addressing post-fire emergency actions when needed. Lastly, remote sensing plays a crucial role in analyzing the evolution of burned areas over time. This can be focused on multiple elements of the ecosystem but is generally focused on soil and vegetation. The assessment of post-fire trajectories is necessary to identify the areas where post-fire recovery is not satisfactory and for the implementation of restoration strategies.

The remote sensing discipline is rapidly advancing thanks to the increasing availability of sensors, data, techniques, and processing capabilities. Thus, in this Special Issue, "Remote Sensing in Forest Fire Monitoring and Post-fire Damage Analysis", we compiled 10 studies [1–10] representing significant advances in the remote sensing of fires, with regard to the different aspects and temporal stages exposed above. In relation to the pre-fire stage, our Special Issue provides new insights into the relevant predictors of fire activity, such as live fuel moisture content [3] and surface fuel load [4], or soil moisture

**Citation:** Fernández-García, V.; Calvo, L.; Suárez-Seoane, S.; Marcos, E. Remote Sensing Advances in Fire Science: From Fire Predictors to Post-Fire Monitoring. *Remote Sens.* **2023**, *15*, 4930. https://doi.org/ 10.3390/rs15204930

Received: 20 September 2023 Accepted: 7 October 2023 Published: 12 October 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

availability [5]. Moreover, Stoyanova et al. [5] explored the potential of land surface temperature status and dynamics as novel indicators of fire risk. In relation to the second stage, a new burned area mapping algorithm using Sentinel-2 [2] is presented, which has been revealed as the most accurate non-commercial imagery for burned area mapping. Furthermore, Park et al. [6] presented a novel approach to improve disaster responses, based on the application of deep learning to detect the multiple elements that might interact in a fire situation. Regarding the assessment of burned areas and fire impacts, our Special Issue includes the first analysis of trends in burn severity at the global scale, revealing the aggravation of fires in many forest biomes [1]. Moreover, Silva Cardoza et al. [7] proposed an improved burn severity algorithm by combining relative phenological correction with former burn severity metrics. In terms of analyzing post-fire trajectories, our Special Issue encompasses a variety of advances, providing cutting-edge information on the drivers and dynamics of post-fire regeneration [8], the performance of physical-based models to measure forest resilience to fire [9], and the identification of optimal parametrizations and wavelengths for LiDAR classifications in post-fire environments [10].

The work provided in this Special Issue contains examples of the multiple advancements in the remote sensing discipline. For instance, the presented studies demonstrate advancements using different remote sensing platforms exploiting imagery from geosynchronous orbit satellites (METEOSAT) [5], sun-synchronous polar orbit satellites (MODIS) [1,3], low-earth-orbit satellites (Sentinel-2) [2,7–9], aircrafts [4,10], or UAVs [6]. Likewise, examples are also provided of how the remote sensing and fire sciences can be advanced using different sensor types (passive [1–3,5–9] and active [4,10]) and methodological approaches (deep learning, machine learning, radiative transfer models, spectral mixture analysis, interpolation techniques, linear models, and spectral indices).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Global Patterns and Dynamics of Burned Area and Burn Severity**

**Víctor Fernández-García 1,2,\* and Esteban Alonso-González <sup>3</sup>**


**Abstract:** It is a widespread assumption that burned area and severity are increasing worldwide due to climate change. This issue has motivated former analysis based on satellite imagery, revealing a decreasing trend in global burned areas. However, few studies have addressed burn severity trends, rarely relating them to climate variables, and none of them at the global scale. Within this context, we characterized the spatiotemporal patterns of burned area and severity by biomes and continents and we analyzed their relationships with climate over 17 years. African flooded and non-flooded grasslands and savannas were the most fire-prone biomes on Earth, whereas taiga and tundra exhibited the highest burn severity. Our temporal analysis updated the evidence of a decreasing trend in the global burned area (−1.50% year−1; *<sup>p</sup>* < 0.01) and revealed increases in the fraction of burned area affected by high severity (0.95% year−1; *p* < 0.05). Likewise, the regions with significant increases in mean burn severity, and burned areas at high severity outnumbered those with significant decreases. Among them, increases in severely burned areas in the temperate broadleaf and mixed forests of South America and tropical moist broadleaf forests of Australia were particularly intense. Although the spatial patterns of burned area and severity are clearly driven by climate, we did not find climate warming to increase burned area and burn severity over time, suggesting other factors as the primary drivers of current shifts in fire regimes at the planetary scale.

**Keywords:** fire severity; burn severity; spatial patterns; trends; biomes; continents; climate warming

#### **1. Introduction**

Fire activity plays a key role in shaping ecological communities, biogeochemical cycles, climate, and human lives and assets [1,2]. More than half of the land surface on Earth is prone to fire [3], with about 3% burning annually [4]. As a result, landscape fires generate annual emissions estimated in 2.2 PgC, an equivalent of 23% of global CO2 from fossil fuels [5], influence global albedo [6], and cause premature deaths from poor air quality, dozens of direct fatalities, and annual economic losses estimated at above US\$2500 million [7]. Despite these global figures, burned area (BA) and its directly related fire regime variable, the fire frequency [8], are largely heterogeneous across space and time because of differences in their main determining factors. Among the BA determining factors there are fuel load, which is linked to primary productivity and herbivory [1,9,10]; fuel connectivity, which is enabled by non-fragmented landscapes and non-sparse vegetation [11,12]; flammability of fuels, closely linked to weather and vegetation properties [2,13,14]; and ignition sources, which might be natural or anthropogenic, the last exceeding 90% of ignition events in most terrestrial biomes [8]. In the same way, the ecological consequences of fire vary depending on the intrinsic ecosystem traits, post-fire environmental conditions, and fire characteristics, mainly burn severity (BS). BS is closely linked to fire intensity and here

**Citation:** Fernández-García, V.; Alonso-González, E. Global Patterns and Dynamics of Burned Area and Burn Severity. *Remote Sens.* **2023**, *15*, 3401. https://doi.org/10.3390/ rs15133401

Academic Editor: Luis A. Ruiz

Received: 31 May 2023 Revised: 27 June 2023 Accepted: 2 July 2023 Published: 4 July 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

is defined as the immediate degree of overall environmental change caused by a fire in ecosystems, including biomass loss, vegetation mortality, and biochemical and physical impacts on soil [15]. Moreover, BS determines the post-fire processes, playing a key role in the future ecosystem pathways [15–18], and its characterization is essential for the refinement of carbon emission models [5,15,19].

Nowadays, the exploitation of some types of satellite imagery provides globally consistent BA products based on the detection of active fires and the spectral changes in surface reflectance [4,20]. Remote sensing analysis revealed grasslands and savannas as the most fire-prone biomes on Earth, mainly those located in Africa and Northern Australia. In these regions many areas burn annually or biennially, exhibiting a large quantitative difference from the rest of the biomes, where generally less than 2% of the area burns every year [4,5]. Satellite imagery also allows the characterization of spatial patterns of BS [19,21,22]. BS has been customarily assessed by spectral indices based on the near-infrared and shortwave infrared reflectance, such as the difference in the normalized burn ratio (dNBR) [16,22,23]. The dNBR accurately matches field measurements of the overall environmental change caused by a fire in multiple ecosystems (mainly calculated through composite indices combining several ecosystem metrics) and has been adopted as a standard BS metric by the EFFIS and MTBS programs of the European Union and the United States of America, respectively [19]. Despite the lack of comprehensive global spatiotemporal BS characterizations [7], current knowledge of maximum fire intensities [20], biomass burned fractions [5] and BS drivers (mainly the fuel load available to burn and fuel continuity) [24–26] anticipate a large spatial variability in BS across the globe.

The evolution of BA and BS arouses great interest in the media and in the scientific community, which frequently warned about the increase or worsening of fire activity during recent years [7], often attributing that assumed trend to climate change (see examples in [17,27–29]). In this sense, climate warming is expected to cause increases in fire weather danger in many regions [14], which is a driver of BA in a large proportion of Earth's land surface [13] and influences BS [24,25]. However, former empirical evidence of BA and BS increases are often based on constrained observations, in terms of timescales or spatial coverage [7]. In fact, global quantitative BA analysis has shown significant decreases in the global BA between 2003 and 2015, mainly concentrated in the tropical savannas of Africa, South America, and Asian steppes, albeit BA increases have been detected in many regions such as the Siberian Arctic [30] and others dominated by closed canopy forests [4], including Australia where BA increases have been fostered by climate change [31]. Likewise, little is known about BS trends at the global scale as most of the literature focuses on USA, Australia, and Mediterranean Europe, which might lead to a global "Western" biased perception [7]. For instance, it is well documented a shift from low to high severity fire regimes in southwestern US forests, caused by the implementation of fire suppression policies after the European settlement [1], and extensive spatiotemporal studies have revealed a generalized increase in high-severity fires in some parts of Western US between 1984–2015 [21]. In Australia, an increase in the proportion of BA at high severity has been detected between 2013 and 2017 [17], and in Europe, an increased prevalence of extreme wildfire events attributed to climate change and human alterations of landscapes has been reported [29].

Here, we characterize the spatial patterns and temporal trends of BA, BS, and BA by BS levels at the global scale and by biomes and continents for the period 2003–2019 (17 years). In addition, we studied the potential relationships between spatiotemporal patterns of BA, BS, and BA by severity levels with climate variables (mean annual temperature, annual temperature range, annual precipitation, and annual precipitation range) to identify the former role of climate change on fire activity. To achieve these objectives, we used NASA-MODIS-derived products at 500 m spatial resolution and ERA5 ECMWF re-analysis climatic data.

#### **2. Materials and Methods**

#### *2.1. Data Sources*

The BA and BS data were obtained from the MOSEV database [19]. MOSEV has been developed based on MODIS information from the MCD64A1 collection 6 [32], which is the standard NASA global BA product and probably the most used by the scientific community [19]; and from the Terra MOD09A1 and Aqua MYD09A1 version 6 products, which provide the highest quality surface reflectance observations in eight-day periods. MOSEV offers, among other information, monthly global data at 500 m spatial resolution from 2000 to 2019 of date of burning and BS measured by the dNBR spectral index ranging from −2000 to 2000 [19]. Burn severity in the MOSEV database refers to short-term BS, also called initial assessment of BS as the closest pre- and post-fire MOD09A1 and MYD09A1 scenes are used for calculations. The alternative to initial BS assessments, are the extended assessments that are often made one year after burning and thus are influenced by ecosystem recovery [15,16], BS in the MOSEV database is directly related to short-term BS calculated with Landsat imagery (R = 0.74 for dNBR; R = 0.42 for RdNBR) [19]. Here, we used the dNBR instead of the relativized index RdNBR, also available in MOSEV, according to (I) its better performance when using MODIS; (II) to its similar or better fit to field composite metrics of BS such as the Composite Burn Index in different types of ecosystems; and (III) to be consistent with our definition of BS, which is understood as the overall degree of environmental change caused by fire [19]. For the present study, the period 2003–2019 was selected for consistency because Aqua products were not available for entire years before 2003. Mean temperature and precipitation data were extracted for each month of the time series from the ERA5 ECMWF reanalysis. ERA5 was downloaded from the Copernicus Climate Data Store (CDS), using the CDS Python API.

#### *2.2. Data Preparation*

All data were extracted globally and by the regions conformed by the intersection of terrestrial biomes (i.e., global vegetation units) [33] with the continent boundaries [34]. The extraction of the global and regional data from both MOSEV and ERA5 was designed in the R programming language [35]. All the routines were implemented in the supercomputing facilities of the Spanish Research Council (CSIC).

For each year, we computed globally and by regions the BA, the percentage of the land burned, the mean BS, the percentage of land burned at different BS levels, and the percentage of the BA burned at different BS levels. The number of BS levels and threshold values depend on the user's objectives [16,19,36]. The differentiation of three categories (low, moderate, and high BS) is probably the most common approach [19,37,38] and many studies have followed the pioneering proposal by Key and Benson [16] to differentiate levels. In the present study, we differentiated low, moderate, and high BS levels by the standard <270, 270–440, and >440 dNBR ranges, respectively, following the thresholds proposed in [16]. We have used the dNBR value of >440 to differentiate the high severity rather than other common higher values (e.g., 660) as those have been found to be excessive by multiple studies [37,38], and because values above 660 are scarce in the MOSEV database [19]. Moreover, we highlight that we have applied the same thresholds to all the regions for consistency; thus, same categories reflect similar overall spectral change with respect to the pre-burn situation regardless of the region.

Likewise, for each year we computed the mean annual temperature, the annual temperature range as the difference between the hottest and coldest months, the annual precipitation, and the annual precipitation range as the difference between the wettest and driest months. This raw database with the annual values of all study variables was checked for coherence before performing the statistics based on regional data. Thus, biomes corresponding to lakes, water bodies, rock, and ice were removed, as well as those regions with BA registered less than 10 years of the study period (<60% of the time).

#### *2.3. Data Analysis*

All the statistical analyses were performed in R software [35]. The spatial patterns of all fire variables (BA, percentage of the land burned, mean BS, percentage of land burned at different BS levels, percentage of the BA burned at different BS levels), and climate variables (mean annual temperature, annual temperature range, annual precipitation, and annual precipitation range) were studied by regions calculating the mean (±standard error) values for the entire study period (2003–2019). Likewise, we computed the mean global values of all fire variables.

To perform the temporal trend analyses (2003–2019; n = 17), we normalized the global and regional annual data of fire variables to the corresponding 2003 values (with values of 2003 set to 100%), to facilitate interpretation and intercomparison. Afterward, we calculated the Sen robust regression estimator with 95% confidence intervals using the zyp.sen and confint.zyp functions of the zyp package [39], and the Mann–Kendall significance using the cor.test function. The Sen robust regression estimator is a non-parametric statistical analysis computed from the median of the slopes of all lines through pairs of points, and thus it is insensitive to outliers. The non-parametric Mann–Kendall significance indicates if a variable consistently decreases or increases over time.

The influence of climate variables (independent) on each fire variable (dependent) was studied with the data by regions (n = 64; n = 63 for BA at high severity analyses because of the absence of this BS level in one region) using two complementary statistical approaches: first, the relationship of each single climate variable with each fire variable was shown by fitting the local polynomial regression (loess), and the Spearman's rank correlation coefficients and significances were calculated using the cor.test function; secondly, conditional inference regression trees were implemented to account for complex combined effects of climate variables. All mean 2003–2019 BA data were square root transformed for these statistical analyses and to facilitate outputs visualization. Regression trees were built with the ctree function of the partykit package [40], including as predictors the four climate variables together. Regression trees include the most significant partitions up to *p* < 0.10 calculated with the approximated finite-sample distribution of Monte Carlo for each node. The variance explained by the regression tree predictions against our data was also shown. The analyses were run to analyze both, the synchronous relationships between fire and climate by using the 2003–2019 means; and the influence of temporal trends, using the Sen slopes as input data.

To analyze the relationships between pairs of fire variables, as well as between pairs of climate variables, correlation matrices with the loess fit and Spearman's rank correlation coefficients and significances were computed for the synchronous spatial data (mean 2003–2019 values) and trend data (Sen slopes) by regions (n = 64) using the pairs.panels function of the psych package [41].

#### **3. Results and Discussion**

#### *3.1. Spatial Patterns of Burned Area and Burn Severity*

We observed that 4.78 ± 0.12 Mkm2 (mean ± standard error), and around 3.26 ± 0.08% of Earth's land surface burned annually in the period 2003–2019 (Table 1), similar to NASA-MODIS MCD64A1 BA data reported for 2003–2015 [4]. The spatial patterns of absolute BA were highly heterogeneous among regions (biomes and continents) (Table A1). The largest contributors to global BA were the tropical and subtropical grasslands, savannas, and shrublands of Africa (2.79 ± 0.07 Mkm<sup>2</sup> year<sup>−</sup>1), Australia (0.36 ± 0.03 Mkm<sup>2</sup> year<sup>−</sup>1) and South America (0.26 ± 0.02 Mkm2 year−1), a consequence of both their proneness to fire and their vast extent.

Relativizing the BA to the total extent of each region (Figure 1; Table A1), we found that the most fire-prone regions were the flooded grasslands and savannas of Africa (26.97 ± 1.01% year−1), as well as the tropical and subtropical grasslands, savannas and shrublands of Africa (20.07 ± 0.49% year<sup>−</sup>1), and Australia (17.01 ± 1.32% year<sup>−</sup>1). There, the seasonally dry climate enables positive feedback interactions between primary productivity, fuel connectivity, and fire [1,12,42] which are essential for maintaining the open structure of these biomes [9,43]. Among grasslands and savannas, the highest potential to burn in those seasonally flooded has been attributed to greater fuel accumulation rates as consequence of a higher productivity and the low herbivory of their less palatable vegetation [10]. We found other biomes to be particularly fire-prone with more than 5% of their area burned annually (Figure 1; Table A1), including the tropical and subtropical dry broadleaf forests of Africa and montane grasslands and shrublands of Africa and Australia. Likewise, a surprisingly high BA fraction (4.38% to 3.42% year<sup>−</sup>1) was detected in tropical and subtropical moist broadleaf forests of Africa and Australia as fire is not an intrinsic ecological process in rainforests [1,31]. The BA fraction was below 1.20% year−<sup>1</sup> in all temperate forests and boreal biomes, except for the Australian temperate broadleaf and mixed forests (1.72% year−1). Although Mediterranean biomes are considered among the most fire-prone on Earth [3], their BA ranged from 0.45% year−<sup>1</sup> in Europe to 1.56% year−<sup>1</sup> in North America. Our results revealed that most biomes burned more in Africa and Australia than on the other continents. In Africa, extensive burning is facilitated by the low population density and decreases in grazing by bulk-feeding herbivores [12], whereas the proneness of Australian biomes to fire is fostered by the extreme climate seasonality, the high number of dry lightning events, aboriginal fire management practices, the physiognomy of eucalypts vegetation [31,42,44], and even by the colonization of exotic grasses [45].

**Table 1.** Global mean values and normalized trends (value in 2003 = 100%) in burned area (BA), burn severity (BS), and BA by BS levels at the global scale for the period 2003–2019. Both, BA, and BA by BS levels are presented in absolute area values and in percentage with respect to the total extent of the global land surface, BA by BS levels are also presented in percentage with respect to the global BA. Note that BA and BA by BS levels, whether calculated as absolute area or as percentage with respect to the total extent of the region results in the same trend values. The entire land area of the world was considered in calculations, including regions with scarce or null fire occurrence. Significant trends (*p* < 0.05) are denoted by boldface and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001). SE: standard error; CL: confidence limit.


BS spatial patterns were closely linked to biomes, being more consistent among continents than BA (Figure 1, Table A1). Generally, the BS patterns detected in the present study are in line with those of fire radiative power emitted by fires detected in previous work [20], revealing at the global scale the assumed relationship between BS and fire intensity [15]. In general, we found the highest mean severities in the taiga, followed by tundra, temperate coniferous forests, and Mediterranean biomes. BS was particularly high in North American taiga, which is characterized by high-intensity stand-replacing crown fire regimes [25], whereas lower mortality (42% lower) and combustion completeness (36% lower) characterize the Eurasian boreal forests, due to different traits of the dominant

species [2,46]. The high severity of tundra fires is a consequence of the burning completeness of vascular vegetation and moss cover, as well as of the upper soil organic strata accumulated for decades, which contributes to permafrost degradation [47]. It is important to note that biomes at high latitudes have the largest topsoil carbon stocks on Earth (0–30 cm) [48], and seasonally become available to severely burn. In Mediterranean biomes, high severity results from a marked crown fire regime [1,3]. The lowest severities were found in tropical coniferous forests, mangroves, and in deserts and xeric shrublands; in the first, combustion completeness is flammability-limited, and the last BS is constrained by fuel-limitations. Tropical and subtropical forests in Africa exhibited lower BS than in Asia and North America, which can be a consequence of burning at smaller fire patches [46] as fire size covariates with BS [17], and potentially because of the predominance of burning in lower tree cover ranges (25–75%) [46]. Low severity values were also found in non-flooded grassy biomes, and intermediate severities generally corresponded to temperate forests. These global patterns confirmed the tradeoffs between BA fractions and BS (ρ = −0.25; *p* < 0.05) (Figure 1) previously detected for fire intensity [20], at the time that supports the fuel load available to burn as the primary driver of BS [25,26,49].

The analysis of BA by different severity levels, particularly at high severity, is useful to acknowledge the territorial magnitude of fire impacts [21]. Globally, we revealed that 0.21 ± 0.00 MKm<sup>2</sup> (0.14 ± 0.00%) of land surface burned annually at high severity (Table 1). In absolute terms, the extent of BA affected by high severity was significantly linked to the BA extent (ρ = 0.87; *p* < 0.01) (Figure 1), the largest contributors to the global high severity BA are the tropical and subtropical grasslands, savannas, and shrublands of Africa, South America, and the Asian taiga (Table A1). However, relativized to the total extent of each region, we found that the regions most affected by high-severity burning were the montane grasslands and shrublands of Australia (2.23 ± 1.77% year<sup>−</sup>1), flooded grasslands and savannas of North America (1.85 ± 0.13% year−1), Africa (1.17 ± 0.12% year−1) and South America (0.96 ± 0.19% year−1), tropical and subtropical grasslands, savannas, and shrublands of Oceania (1.60 ± 1.59% year−1) and South America (0.55 ± 0.04% year−1), and Mediterranean North America (0.52 ± 0.11% year−1) (Figure 1; Table A1). In the rest of the regions, the high severity BA constituted less than 1% of their extent. This is a direct consequence of the BA in each region (ρ = 0.71; *p* < 0.01), being related to a lesser extent to the mean BS (Figure 1). Relativizing the BA at high severity to the total BA in each region, similar patterns to BS were detected (ρ = 0.82; *p* < 0.01) (Figure 1), with around 40% of BA exhibiting a high severity in all boreal forests and non-European tundra (Figure 1; Table A1).

#### *3.2. Temporal Trends in Burned Area and Burn Severity*

World BA decreased −1.50% year−1; *<sup>p</sup>* < 0.01 which in absolute terms is equivalent to a decrease of 1.22 MKm<sup>2</sup> between 2003 and 2019 (25.5% less BA with respect to 2003) (Table 1). These results elucidate an intensification of the decline in BA compared to the results already reported for the period 2003–2015 [4]. Analyzing BA trends by regions (Figure 2; Table 2), we found significant BA decreases in the European taiga (−6.21% year<sup>−</sup>1; *<sup>p</sup>* = 0.04), temperate grasslands, savannas and shrublands of Asia (−2.77% year<sup>−</sup>1; *<sup>p</sup>* = 0.01), Mediterranean Europe (−2.31% year−1; *<sup>p</sup>* = 0.03), tropical and subtropical dry broadleaf forests of North America (−2.11% year−1; *<sup>p</sup>* < 0.01), montane grasslands and shrublands of Africa (−2.08% year−1; *<sup>p</sup>* < 0.01) and Asia (−1.93% year−1; *<sup>p</sup>* = 0.02), and in most of tropical Africa including grasslands, savannas and shrublands (−1.64% year−1; *<sup>p</sup>* < 0.01), moist broadleaf forests (−1.57% year<sup>−</sup>1; *<sup>p</sup>* = 0.02) and mangroves (−1.37% year<sup>−</sup>1; *<sup>p</sup>* = 0.02). Generalized decreases in BA have been formerly attributed to a decrease in the number of ignitions and to a lesser extent, to decreases in fire size, driven by human activity [4]. In this sense, the increase in human population, cropland area, and livestock density cause decreases in fire activity in the fire-prone open biomes [4], whereas increased efficiency in fire prevention, detection, and extinction, and abandonment of fire use in agriculture contribute to the decreasing trend in other regions [7,50].


**Figure 1.** Mean values of burned area (BA), burn severity (BS), and BA by BS levels by regions (biomes × continents) for the period 2003–2019. BA values are represented by square root scaled bars and expressed in percentage with respect to the total extent of the regional land surface. The percentage of the BA burned at each BS level is proportional to their fraction within the BA bars. BS is represented by red points and expressed in dNBR units ranging from −2000 to 2000. S.Am: South America, Oce: Oceania, N.Am: North America, Eu: Europe, Aus: Australia, As: Asia, Afr: Africa.

**Figure 2.** Normalized trends (value in 2003 = 100%) in burned area (BA) by regions (biomes × continents) for the period 2003–2019. The color ramps are square root scaled.

Globally, BS exhibited a non-significant increase (0.13% year−1; *p* = 0.34) (Table 1). However, we found a large continental disparity, and regions with significant increases outnumbered those with significant decreases in BS (Figure 3, Table 2). The largest increases were found in South America, including temperate grasslands, savannas, and shrublands (2.47% year−1; *p* = 0.03) and all forest biomes (1.49 to 3.56% year−1; *p* < 0.05) except mangroves. Likewise, we detected significant increases in BS in tropical and subtropical moist forests of Australia (2.23% year−1; *p* = 0.04), Mediterranean Africa (2.16% year−1; *p* < 0.01), and Asian taiga (2.13% year−1; *p* = 0.03). Land cover changes might explain BS trends in many regions, for instance, most of Mediterranean Africa, Central and Southern Chile, and Siberian taiga have experienced long-term woody encroachment and forest expansion [51], a fuel accumulation that in these regions lead to more severe fires [3,49]. Likewise, BS increases in the tropical rainforest of South America may respond to several causes. First, to increases in fuel continuity in the Eastern Brazilian coast, which is frequently burned [50]; second, to the loss of primary forests [51] that are characterized by no fire or low severity surface fires [1]; and third, to decreases in lower atmosphere moisture boosted by these forest losses that resulted in increased droughts the last years of the study period [52], making large fuel loads available to burn. Significant decreases in BS were only found in montane grasslands and shrublands of Asia (2.94% year−1; *p* = 0.04) and in the European taiga (2.56% year<sup>−</sup>1; *p* = 0.04), which inversely to the Asian, has experienced decreasing trends in fuel continuity and increasing moisture which have detrimental effects on BA [50] and likely on BS [24,26].

The analysis of BA by BS levels showed that decreases in global BA were mainly at the expense of decreasing low (−1.62% year−1; *<sup>p</sup>* < 0.01) and moderate severity BA (−1.27% year<sup>−</sup>1; *<sup>p</sup>* < 0.01) (Table 1). Globally, we found trends in BA at high severity (−0.63% year<sup>−</sup>1; *<sup>p</sup>* = 0.09) to be positively related to trends in BA (<sup>ρ</sup> = 0.67; *<sup>p</sup>* < 0.01) as well as to trends in BS (ρ = 0.40; *p* < 0.01) (Figure 2). Regionally, we detected the highest increases in BA at high severity in temperate broadleaf forests of South America (40.20% year−1; *p* = 0.04) (Figure 4A; Table 2), which can be a consequence of the expansion of exotic *Pinus* and *Eucalyptus* plantations with associated fuel load increases up to 40 Mg ha−<sup>1</sup> between 1999 and 2006, accompanied by a large drought-driven intensification of fire activity between 2010 and 2015 [53]. We observed the BA at high severity to also aggravate in Australian tropical and subtropical moist forests (28.02% year−1; *p* < 0.01), which are vulnerable to fire, as it favors alternative open biome states [44]. BA at high severity also increased in Australian grasslands, savannas, and shrublands (6.85% year−1; *p* = 0.01). In Australia, increases in fuel continuity [50] and in fire weather [31] in the last decades have been detected, which, along with the unprecedented fire events that occurred in 2019, might contribute to the detected patterns. We found significant increases in BA at high severity in the deserts and xeric shrublands and tropical and subtropical dry broadleaf forests of Asia that locally can be attributed to encroachment [49] and to increases in fuel continuity [50]. The largest decreases in high severity BA were observed in European taiga (−6.03% year−1; *<sup>p</sup>* = 0.01) and in tropical and subtropical grasslands, savannas, and shrublands of North America (−5.96% year<sup>−</sup>1; *<sup>p</sup>* = 0.04), a consequence of decreases in both BA and BS. Moreover, a significant global increase in the fraction of BA that is burned at high severity was detected (0.95% year−1; *p* = 0.04), and regionally, trends in the fraction of BA burned at high severity were closer to BS (ρ = 0.51; *p* < 0.01) than to BA (ρ = 0.11; *p* = 0.38) (Figure 2), confirming the aggravation of impacts within burned areas in a vast proportion of South America, and in large parts of Northern Australia and Asia (Figure 4B, Table 2). In the tropical forests of Australia and the Amazon, it has been revealed that contemporary logging regimes and silvicultural practices exacerbate burn severity [54,55], suggesting decision makers and forest managers have a determinant role in past, present, and future fire regimes.

**Figure 3.** Normalized trends (value in 2003 = 100%) burn severity (BS) by regions (biomes × continents) for the period 2003–2019. The color ramps are square root scaled.

#### *3.3. Relationships with Climate Variables*

The climate spatial patterns (Figure 3) determined the spatial patterns of fire activity (Figures 5 and 6). Local regressions and correlation analysis (Figure 5) showed monotonic increases in BA (ρ = 0.58; *p* < 0.01), and monotonic decreases in BS (ρ = −0.35; *p* < 0.01) and in the BA at high severity (ρ = −0.54; *p* < 0.01) towards the warmest regions (tropics), in agreement with the results shown by the regression trees (Figure 6A–C). Moreover, the regression trees revealed mean annual temperature as the primary climate driver of regional differences in BA, BS, and BA at high severity (Figure 6A–C). Local regressions (Figure 5) and correlation analysis also showed an opposite pattern for annual temperature range, but regression trees (Figure 6A–C) suggest lower importance of annual temperature range and the other climate variables (*p* > 0.05), probably because of their strong interdependencies (Figure 4).

We found climate trends to be weakly related to all BA and BS trends (Figures 6D and 7), the detected trends being opposite to the assumption of increases in BS with climate warming. Thereby, we detected a significant inverse relationship between climate warming and the proportion of BA at high severity (ρ = −0.26; *p* = 0.04) (Figure 7), and both the local regression and the regression tree (Figure 6D) showed that the most stable regions in terms

of proportion of BA burned at high severity were those experiencing the highest warming (>0.03 K year<sup>−</sup>1). In this context, the long-term outcome of increased temperatures for fire regimes has been disclosed complex as it is caused by shifts in fire weather [31], but also by the warming effects on fuels through changes in productivity, vegetation composition [2,26], and fuel accumulation in the soil uppermost layers that are controlled by productivity– decomposition equilibriums varying at fine scales [56]. Moreover, the climate shifts would have different impacts on fire regimes depending on the former climate, productivity, and induced vegetation trajectories [24,27]. This complexity is expected to cause differential changes in fire regimes across regions, which some authors claim to be scarcely plausible until the coming decades [13,14]. Likewise, these premises highlight the difficulty of linking climate change to generalized increases in BA and BS worldwide, as different variations are expected depending on regional intrinsic characteristics.

**Figure 4.** Normalized trends (value in 2003 = 100%) in burned area (BA) at high severity by regions (biomes × continents) for the period 2003–2019. The map in the upper part (**A**) shows the trends of BA at high severity calculated as absolute burned area at high severity (when calculated as percentage with respect to the total extent of the region the result is the same). The map in the bottom (**B**) shows the trends of BA at high severity calculated as a percentage of the burned area that burned at high severity. The color ramps are square root scaled.

**Figure 5.** Scatterplots showing the synchronous relationships between the climate and fire variables for the period 2003–2019. Red lines show the local polynomial regression fitting (loess), and shaded areas ± 95% confidence intervals. Each panel also shows the Spearman's correlation coefficient (ρ), and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001). Note that the BA axes are square root transformed. BA: burned area, BS: burn severity, T: mean annual temperature: T range: annual temperature range, P: annual precipitation, P range: and annual precipitation range.

**Figure 6.** Conditional inference regression trees showing the most important climatic variables in determining the spatial patterns of burned area (BA) (**A**), burn severity (BS) (**B**), BA at high severity (**C**), and normalized trend (2003 = 100%) in the percentage of the BA burned at high (**D**). No further trees were able to grow for the other studied fire variables and trends at the selected confidence level for recursive binary partitioning (*p* < 0.10). R<sup>2</sup> values on the bottom of panels represent the variance explained by the regression tree predictions on our dataset. Red dots in the boxplots indicate the mean values. T: mean annual temperature, T range: annual temperature range, P: annual precipitation.

**Figure 7.** Scatterplots showing the relationship between the studied climate trends and normalized trends (value in 2003 = 100%) in fire variables for the period 2003–2019. Red lines show the local polynomial regression fitting (loess), and shaded areas ±95% confidence intervals. Each panel also shows the Spearman's correlation coefficient (ρ), and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001). BA: burned area, BS: burn severity, T: mean annual temperature: T range: annual temperature range, P: annual precipitation, P range: and annual precipitation range.

#### *3.4. Implications and Final Considerations*

Our study updates BA data and constitutes the first global assessment of BS spatiotemporal patterns. The information, provided by biomes and continents, is essential in revealing increasing trends in BS in several regions and increases in the fraction of BA that is burned at high severity at the planetary scale, providing scientific support for widespread assumptions. However, we did not find evidence in line with the hypothesis that climate change is increasing the mean BS worldwide. Nevertheless, several issues should be considered when interpreting our findings. First, the used BA and BS data at 500 m spatial resolution are appropriate for global fire analysis and consistent for studying

trends [4,19,31], but have limitations in terms of quantifying BA absolute values. Accordingly, it has been demonstrated that the use of finer spatial resolution imagery at regional scales can result in larger BA, particularly in those regions with a predominance of small fires [57,58]. For this reason, as sensors and computational capabilities improve, we encourage future work to analyze BA and BS trends using higher spatial resolution imagery or statistically refined BA data [58], and even different methods to quantify burn severity [59]. Secondly, although there are biome–fire regime relationships, it is important to note that most biome regions have been largely modified by human activities such as agriculture, livestock, or urbanization, which also have an important role on fire regimes [2,4,8,20,60]. In this sense, we recommend exploring the implications of these drivers in determining the results presented in this study. Third, our study period is limited to 17 years because of no prior availability of both MODIS Aqua and Terra sensors. This could favor interferences of decadal ocean and atmospheric oscillations [61] in our results, but we assume little influence of El Niño Southern Oscillation (ENSO) due to 2–7 years oscillation of this phenomenon and the balanced events along our study period. Fourth, despite not finding strong evidence of climate change aggravation of BA and BS, we confirmed that climate is an essential variable influencing fire regimes worldwide [1,13,27], and further studies at finer scale map units of analysis are advisable to further disentangle the complex climate–fire–human interactions in the current context of climate change.

#### **4. Conclusions**

This study used 17 years of MODIS-derived data to analyze the spatiotemporal patterns of BA and constitutes the first BS analysis at the global scale. In addition, we showed the spatiotemporal patterns of BA and BS by regions (biomes and continents), relating them to several climate variables.

Our results updated the already known spatial patterns of BA across the globe and corroborated significant decreases in global BA. Our triple-way analyses of BS trends detected that globally (I) the mean BS exhibits non-significant increases; (II) the fraction of land affected by high BS shows non-significant decreases, as it is linked to the decreases in BA; and (III) the fraction of burned area that is affected by high BS is significantly increasing. In addition, the number of regions showing significant increases in mean BS and burned area at high BS outnumbered those with significant decreases.

We found close relationships between the spatial patterns of fire variables (BA and BS metrics) and climate variables (temperature, precipitation, and their respective interannual ranges), but our analysis by regions has not found evidence that climate warming is increasing BA nor BS, suggesting other factors as the primary drivers of change.

Given the great technical and computational advances that are currently taking place, future work may continue our analysis by using higher-resolution images, smaller analytical units, longer periods, or different methods of quantifying BS. We also encourage future research to analyze region-by-region the implications that our results may have in the field of ecology, climate regulation, and in the effectiveness and design of environmental management strategies.

**Author Contributions:** Conceptualization, V.F.-G. and E.A.-G.; methodology, V.F.-G. and E.A.-G.; formal analysis, V.F.-G. and E.A.-G.; investigation, V.F.-G. and E.A.-G.; writing—original draft preparation, V.F.-G.; writing—review and editing, V.F.-G. and E.A.-G.; visualization, V.F.-G. and E.A.-G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present study has been conducted without project funding. Víctor Fernández-García is supported by a Margarita Salas post-doctoral fellowship from the Ministry of Universities of Spain, financed with European Union-NextGenerationEU funds and granted by the University of León. Esteban Alonso-González has been funded by the Centre National d'Etudes Spatiales (CNES) postdoctoral fellowship.

**Data Availability Statement:** Data will be made available upon request to the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.




**Table A1.** *Cont.*


**Table A1.** *Cont.*

**Figure 1.** Synchronous relationships between the spatial patterns of mean annual burned area (BA) expressed in mean absolute values (Mkm2 yr−1) and in percentage with respect to the total extent of each region (% land yr−1), burn severity (BS) expressed in dNBR units ranging from <sup>−</sup>2000 to 2000, and BA at high severity expressed in absolute values (Mha yr<sup>−</sup>1), in percentage with respect to the total extent of the land extent in each region (% land yr−1), and in percentage with respect the total BA in each region (% BA yr<sup>−</sup>1) for the period 2003–2019. Red lines in the scatterplots show the local polynomial regression fitting (loess) and shaded areas ±95% confidence intervals. Values in the panels are Spearman's correlation coefficients (ρ) between pairs of variables, and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001). Note that all axis for BA variables were square root scaled.

**Table 2.** Mean trends in annual burned area (BA), in burn severity (BS), and in high severity BA (dNBR > 440) by regions (biomes × continents) for the period 2003–2019. Note that trends are normalized (value in 2003 = 100%) and expressed in percentage change per year, and therefore both, trends in mean absolute values of BA (originally measured in Mha year<sup>−</sup>1), and in BA with respect to the total extent of global land extent (originally measured in % land year<sup>−</sup>1) are the same. Slopes were calculated using the Sen slope estimator and significances with the Mann–Kendall test. Significant decreases (*p* < 0.05) are denoted by green boldface, significant increases by red boldface, and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001). Trop.: tropical and subtropical. B: broadleaf forest. C: coniferous forest. M: mixed forest. G.: grasslands. S: savannas. Sh.: shrublands.


#### **Table 2.** *Cont.*



**Table 2.** *Cont.*

**Figure 2.** Relationships between the normalized trends (value in 2003 = 100%) in annual burned area (BA), burn severity (BS), and BA at high severity. Note that BA and BA by BS levels whether calculated as absolute area BA or as percentage with respect to the total extent of the region lead to same trend values. Red lines in the scatterplots show the local polynomial regression fitting (loess) fit and shaded areas ±95% confidence intervals. Values in the panels are Spearman's correlation coefficients (ρ) between pairs of variables, and the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001).

**Figure 3.** Spatial patterns (**A**–**D**) and trends (**E**–**H**) in mean annual temperature, annual temperature range, annual precipitation, and annual precipitation range for the period 2003–2019. Trends were calculated using the Sen estimator and significances with the Mann–Kendall test.

**Figure 4.** Synchronous relationships between the spatial patterns of mean annual temperature (T), annual temperature range (T range), annual precipitation (P) and annual precipitation range (P range) for the period 2003-2019. Red lines in the scatterplots show the local polynomial regression fitting (loess) and shaded areas ±95% confidence intervals. Values in the panels are the Spearman's correlation coefficients (ρ) between pairs of variables. the number of asterisks denotes the level of statistical significance (*p* < 0.05, *p* < 0.01, and *p* < 0.001).

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Post-Fire Forest Vegetation State Monitoring through Satellite Remote Sensing and In Situ Data**

**Daniela Avetisyan 1,\*, Emiliya Velizarova <sup>2</sup> and Lachezar Filchev <sup>1</sup>**


**Abstract:** Wildfires have significant environmental and socio-economic impacts, affecting ecosystems and people worldwide. Over the coming decades, it is expected that the intensity and impact of wildfires will grow depending on the variability of climate parameters. Although Bulgaria is not situated within the geographical borders of the Mediterranean region, which is one of the most vulnerable regions to the impacts of temperature extremes, the climate is strongly influenced by it. Forests are amongst the most vulnerable ecosystems affected by wildfires. They are insufficiently adapted to fire, and the monitoring of fire impacts and post-fire recovery processes is of utmost importance for suggesting actions to mitigate the risk and impact of that catastrophic event. This paper investigated the forest vegetation recovery process after a wildfire in the Ardino region, southeast Bulgaria from the period between 2016 and 2021. The study aimed to present a monitoring approach for the estimation of the post-fire vegetation state with an emphasis on fire-affected territory mapping, evaluation of vegetation damage, fire and burn severity estimation, and assessment of their influence on vegetation recovery. The study used satellite remotely sensed imagery and respective indices of greenness, moisture, and fire severity from Sentinel-2. It utilized the potential of the landscape approach in monitoring processes occurring in fire-affected forest ecosystems. Ancillary data about pre-fire vegetation state and slope inclinations were used to supplement our analysis for a better understanding of the fire regime and post-fire vegetation damages. Slope aspects were used to estimate and compare their impact on the ecosystems' post-fire recovery capacity. Soil data were involved in the interpretation of the results.

**Keywords:** fire impact; post-fire forest recovery; forest landscapes; vegetation indices; orthogonal transformation; Sentinel-2

#### **1. Introduction**

Forest disturbance cycles are associated with exacerbating responses to climate change [1]. Forest fires have been more frequent and severe in recent decades, especially in areas that have experienced climate change pressures for an extended period of time [2]. Due to climate change, high-temperature anomalies continue to occur, which leads to frequent forest fires [3]. The International Panel on Climate Change (IPCC) puts the Mediterranean and its adjacent lands as amongst the most vulnerable regions to the effects of global warming worldwide [4]. The models issued by IPCC agreed on a clear trend of the thermal regime based on a scenario from 1980–2000. An increase in average surface temperatures, ranging between 2.2 ◦C and 5.1 ◦C, for the period 2080–2100 was prognosed. For the same period, the models indicated pronounced rainfall regime changes showing that precipitation over lands might decrease by about 4% to 27%. Studies performed by the Department of Meteorology at the National Institute of Meteorology and Hydrology at the Bulgarian Academy of Sciences (NIMH-BAS) predicted an increase in the annual air temperature by more than 1.8 ◦C for the coming decades in Bulgaria. This fact increases

**Citation:** Avetisyan, D.; Velizarova, E.; Filchev, L. Post-Fire Forest Vegetation State Monitoring through Satellite Remote Sensing and In Situ Data. *Remote Sens.* **2022**, *14*, 6266. https://doi.org/10.3390/rs14246266

Academic Editors: Elena Marcos, Leonor Calvo, Susana Suarez-Seoane and Víctor Fernández-García

Received: 21 October 2022 Accepted: 7 December 2022 Published: 10 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the risk of forest fire frequency and intensity [5]. During the last two decades in Bulgaria, as well as in other countries of the Mediterranean region, many wildfires occurred and had a significant economic, political, social, and ecological impact [6–9].

In recent years, high resolution (HR) and very high resolution (VHR) optical remote sensing has become widespread concerning monitoring needs, and these strategies provide affordable multitemporal and multispectral pictures of the considered phenomena at different scales. Satellite sensors allow measurement of the impact of fires by comparing pre- and post-fire information. Applications of remote sensing technology related to fire ecology, including fire risk mapping, fuel mapping, burn severity assessment, and post-fire vegetation recovery, are widely discussed and accepted [10–15]. These technologies provide a low-cost, multi-temporal means for conducting local, regional, and global-scale fire ecology research. Moreover, the development of new technologies and techniques resulted in their rapid evolution, thus increasing the accuracy and efficiency of earth observation studies and applications [16,17]. Space and airborne sensors have been used to map burned areas, quantify the impact of fire on vegetation over large areas, and characterize post-fire ecological effects [18–20]. Emphasis has been given to the roles of multispectral sensors, lidar, and emerging Unmanned Aerial System (UAS) technologies [14]. Depending on the purpose of post-fire vegetation recovery observation and study, the assessment is performed based on groups of methods, such as image classification, vegetation indices (VIs), spectral mixture analysis (SMA) [21,22], etc. Remote sensing imagery offers an opportunity for obtaining land use and land cover information through image interpretation and classification. Spectral responses are used in image classification to identify the healthy vegetation in individual pixels [14]. Spectra-based classification approaches are conceptually simple and easy to implement [23]. The type or condition of surface features and their dynamics can be assessed by multi-temporal imaging. This type of analysis is fundamental in remote sensing and is typically called change detection [24]. VIs are the most commonly used method for assessments of vegetation state, including vegetation recovery after natural or anthropogenic disturbances [17,25,26].

Post-fire-related conditions are important for forest vegetation recovery. In this context, mapping of the burned area, representing the burn severity, is a standard technique for monitoring the post-fire effects and forest recovery patterns [17,27–29]. It was found that the differenced Normalized Burn Ratio (dNBR) and its relative form (relative differenced Normalized Burn Ratio (RdNBR)) derived from Landsat data correlate with field measurements of burn severity [30]. NBR is also used for monitoring post-fire regeneration over burned areas in ecosystems. Results showed that as vegetation regenerates, the differences between the burns and the reference area for the vegetation index decrease with time [31]. Detailed studies showed that the NIR-based vegetation indices are most appropriate for accurately assessing vegetation recovery [32]. The NDVI is found to be the most used for post-fire recovery studies, as it could be calculated alone without additional field data collection [26,33,34]. However, due to reaching saturation levels before the point where an ecosystem fully recovers its maximum biomass after disturbance, the forest recovery rate could be overestimated when using NDVI [35,36]. For estimations of variations in chlorophyll content and its changes in vegetation after a fire event, the Modified Chlorophyll Absorption Ratio Index (MCARI2) is used [37]. A quantitative analysis of forest degradation resulting from forest fires is performed by introducing the Normalized Differential Greenness Index (NDGI) [38]. The remotely sensed Moisture Stress Index (MSI) is used for canopy stress analysis and is suitable for monitoring coniferous forests and assessing specific damages that cannot be detected using NIR/R vegetation indices [39]. Spectral indices are also used to estimate other ecological parameters related to vegetation recovery. Such parameters are the Leaf Area Index (LAI) [40], the Forest Recovery Index (FRI), and Fractional Vegetation Cover (FVC) [29].

However, differences in fire severity provoke contrasting plant cover and floristic composition when ecosystems recover after forest fires. A multitude of factors such as climate, initial plant mortality, soil characteristics, the topography of the region, and vegetation composition determine the rate of recovery [41,42]. Moreover, vegetation response to fire and post-fire recovery processes differ in the various biogeographical regions and depend on vegetation type and pre-fire vegetation state [43]. For that reason, post-fire vegetation recovery is a complicated process that cannot be assessed by the application of a single, unified spectral index. Despite the advantages of spectral indices for monitoring post-fire vegetation recovery, there is still no single spectral index suitable for assessment of post-fire disturbances or vegetation recovery processes in every ecosystem, scale, and time-lag condition [44]. Even the NBR and its derivative indices developed specifically for fire-affected areas show varying accuracy under different conditions. The NBR and dNBR are considered advantageous for immediate post-fire monitoring, but their accuracy decreases with increasing temporal distance from the fire event and with the progress of vegetation recovery [45]. On the other hand, the Disturbance Index (DI) [46] is effective in monitoring forest disturbances of different origins and their temporal dynamics [47,48]. Due to the involvement of a larger range of spectral information, the DI is considered more accurate in assessing the recovery of undergrowth in forest ecosystems compared to standard monitoring methods using VIs [19,47]. The DI is based on a linear orthogonal transformation of multispectral satellite images [49,50], which increases its ability to differentiate the three main components: soil, vegetation, and moisture [51]. As a result of a fire, these three components are altered to the greatest extent. The DI more precisely separates the unvegetated spectral signatures closely linked to the stand-replacing disturbance from all other forest signatures [46]. This feature makes the DI particularly suitable for monitoring the dynamics of post-fire vegetation recovery.

This paper dissects the forest vegetation recovery stages during a period of six years after a wildfire in the Ardino region in Bulgaria. The study aims to present the potential for exploiting remotely sensed imagery and respective indices of greenness, moisture, and fire severity from Sentinel-2 to support post-fire observation and forest management with an emphasis on fire-affected territory mapping, vegetation damage assessment, fire and burn severity assessment, and their influence on the ability of vegetation to recover. The study used two groups of spectral indices for monitoring the area affected by the wildfire. The first group encompassed VIs using individual spectral bands for their calculation, and the second group included indices utilizing a larger range of spectral information through the orthogonalization of multispectral data. The Normalized Difference Vegetation Index (NDVI), the Modified Chlorophyll Absorption Ratio Index (MCARI2), and the Moisture Stress Index (MSI) are the indices that belong to the first group, and the Normalized Differential Greenness Index (NDGI), the Normalized Differential Wetness Index (NDWNI), and the Disturbance Index (DI) are the indices from the second group. The study took into account the landscape characteristics of the area influencing the processes occurring in fire-affected forest ecosystems and their post-fire recovery dynamics. Ancillary data about pre-fire vegetation state and slope inclinations were used to supplement our analysis for a better understanding of the fire regime and post-fire vegetation damages. Slope aspects were used to estimate and compare their impact on the ecosystems' resilience, vulnerability, and post-fire recovery capacity. Soil data were involved in the interpretation of the results.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area was situated in the southeastern part of the Rhodope Mountains, near Ardino town. The X and Y coordinates of the centroid were calculated to be 25◦6 7"E and 41◦34 30"N. A significant fire took place on 29 July 2016 in the study area (Figure 1). 2016 was the year with the highest wind speed during the summer months of the study period (2016–2021). The average wind speed in the summer of 2016 ranged between 4.7 m/s (in September) and 6.1 m/s (in August), and the average maximum wind speed was between 7.4 m/s (in September) and 8.1 m/s (in July) [52]. Generally, based on hourly weather simulations over the past 30 years in the study area, the maximum wind speed was observed in March (reaching up to 20 m/s average value) and the minimum windspeed was in September (starting from 3 m/s average value). Days with wind speeds above 12 m/s predominated between February and September, and days with wind speeds above 5 m/s predominated between October and January. The wind direction was mainly from the north and northwest [53].

**Figure 1.** Location of the studied area and land cover change in the pre-fire 2013 (**a**) and post-fire 2016 (**b**), 2018 (**c**), and 2020 (**d**) years. Source: Google Earth Pro—Airbus and Maxar Technologies images.

The area affected by the forest fire was 100 ha. The main tree species were Scots pine (*Pinus sylvestris* L.) and Black pine (*Pinus nigra* Arn.). The endemic vegetation in the region refers to the Thracian province of the European deciduous forest area with the main tree species being *Quercus frainetto and Q.cerris*. On karst terrains, *Q. frainetto*, *Q. cerris,* and *Q.pubescens*, mixed with *Carpinus orientalis, Fraxinus ornus, Syringa vulgaris, Cotinus coggygria,* and *Ostrya carpinifolia* can be found [54]. However, because of erosion processes in the 1950s and the expansion of bare lands, massive afforestation with coniferous species has been performed.

Lithogenic diversity was represented by pre-Paleozoic and Paleozoic metamorphic rocks and phyllitoids covered by a Paleogene volcanogenic–sedimentary complex [54]. The terrain was hilly with steep slopes (>15◦). It influenced the fire regime and is also considered a condition for the development of water erosion. Approximately 70% of the area contained slopes above 15◦. The mean value of the slopes was 16◦ and the maximum was 27◦. The slopes predominantly had east, southeast, and south exposures and define warm and dry conditions for vegetation development. The soils were shallow Lithosols mixed with Rendzinas [54].

The study area has a Continental Mediterranean climate with hot summers and mild winters. The minimum amount of precipitation is in summer, and the maximum amount is in winter. Over the last 40 years, a clear trend towards an increase in both average air temperature and precipitation sum has been observed. The average air temperature change showed a linear trend with an increase of 2.2 ◦C, starting from 10.2 ◦C in 1979 and reaching up to 12.4 ◦C in 2021. Summer is getting hotter, with temperatures in July and August consistently higher for the past 16 years. On the other hand, the precipitation sums in the summer months decreased, especially in August, showing a persistent negative tendency during the last 14 years. The winter is getting warmer too. February has been distinguished by sustained higher average air temperatures since 2012. The precipitation sum has also increased during the winters. Overall, there has been an increase (by 126 mm) in precipitation sum over the past 40 years. This increase is primarily due to increased winter precipitations [53]. The indicated climate changes have led to the transition of the studied area from a Warm-summer Mediterranean climate (Csb), according to the Köppen climate classification, to a Hot-summer Mediterranean climate (Csa). In addition, the slopes in the area, which predominantly had east, southeast, and south exposures, determined warm and dry conditions for the development of vegetation. The observed trends significantly increase the risk of fires in the area, loss of biological diversity, and degradation of ecosystems.

#### *2.2. Characteristics of Climatic Anomalies Observed during the Period of 2016–2021*

Figure 2 shows the mean air temperature and precipitation sum anomalies for the period of 2016–2021 on an annual basis and for the summer months (July, August, and September). As a reference, the period between 1981 and 2010 was used.

**Figure 2.** Mean air temperature (**a**) and precipitation sum (**b**) anomalies for the period of 2016–2021 on an annual basis and for the summer months (July, August, and September).

In terms of mean air temperature, 2019 was the year with the highest positive anomaly (+1.4 ◦C). Amongst the summer months, August was with positive anomalies only. The highest anomalous value (+1.6 ◦C) was recorded in 2017. The highest anomalous value for the summers between 2016 and 2021 was recorded in September 2020 (+2.6 ◦C). Three of the years during the period of 2016–2021 were characterized by positive anomalies in all three summer months. Overall for the three summer months, 2020 was with the highest positive anomaly (+4.7◦C). The lowest positive anomaly, overall for the three summer months, was recorded in 2018 (+0.3◦C). The mean air temperature in July 2018 was 1 ◦C lower in comparison with the reference period (Figure 2a) [52].

Higher sums of annual precipitation were observed during the period of 2016–2021. The wettest year was 2021, which had a 36% higher precipitation sum in comparison with the reference period. There was no definite trend in the distribution of precipitation during the summer months. All three summer months of 2016 and 2020 had lower precipitation. In 2020, the total precipitation for the summer was 106% lower in comparison with the reference period. In September 2020 only, the precipitation sum was lower by 83%. The wettest summer was in 2019. In that summer, the total amount of precipitations was 366% higher in comparison with the reference period (Figure 2b) [52].

#### *2.3. Data*

#### 2.3.1. In Situ Data

In situ data included climatic data, soil data, and field observations.

The climatic data were based on measurements at the Kardzhali meteorological station (331 m.a.s.l) situated 21 km from Ardino [53] and on hourly weather simulations with 30 km spatial resolution over the past 30 years [52]. The in situ climatic data included mean air temperature and precipitation sum for a period between 1979 and 2021 as well as data about the mean air temperature and precipitation sum anomalies for the period of 2016– 2021 on an annual basis and for the summer months (July, August, and September) [53]. In addition, ERA5 model data, which combined satellite and in situ historical observation, were used to outline how climate change has already affected the Ardino region in the last 40 years [52].

The soil data included soil types [54], soil chemical composition, and organic matter content, which were measured in a pine-dominated mixed woodland (790 m.a.s.l) situated near Ardino in 2015 [55]. Soil field data is part of LUCAS 2015 Topsoil datasets, which is freely available through the European Soil Data Centre (ESDAC) of the Joint Research Centre [55]. Because the soil data refer to the pre-fire period and the soil samples were not taken from the affected area specifically, they were used only for result interpretation as auxiliary data.

Interactive three-dimensional panoramas were used as a means of field observations. They were acquired via Google Street View technology in the summer of 2021. These were used to generate high-quality photographs of eight locations affected by the fire [56].

#### 2.3.2. Satellite Data

Satellite data acquired from Sentinel-2A and Sentinel-2B multispectral sensors of the European Space Agency Program for Earth Observation "Copernicus" [57] were used to assess post-fire vegetation recovery. The temporal resolution of every individual Sentinel-2 satellite is ten days, and their combined resolution is five days. More detailed information about the spectral and spatial resolution of the Sentinel-2 satellites can be found in Table A1 [57].

The Sentinel-2 image acquired on 10 July 2016 was representative for the period before the fire event (29 July 2016), and the images acquired on 24 August 2017 and 2018, 14 August 2019, 02 September 2020, and 23 August 2021 were used for assessment of the forest vegetation state after the fire.

High-resolution forest layers (HRLs), which are freely available through the Copernicus Land Monitoring Service [58], were used in the validation process. The layers included

Tree Cover Density (TCD) and Forest Type (FTY) products. The TCD product represents the level of tree cover density in a range from 0–100% for 2012, 2015, and 2018 reference years. The Forest Type product represents the dominant leaf type with a Minimum Map Unit (MMU) of 0.5 ha. Both products are pixel-based, and the minimum mapping width is 20 m. The forests' HRLs in 2015 were verified for the territory of Bulgaria using in situ data [59]. The verification procedure was performed on three levels, including the highly recommended quantitative verification. According to the results obtained by Tepeliev et al. [59], the HRLs used in the present study are generally correctly mapped for the territory of Bulgaria. Hence, we assumed that they could be used as independent reference data in the validation procedure.

A Digital Elevation Model (DEM) with a spatial resolution of 25 m was used to obtain slopes and slopes' aspects. The dataset is freely available through the Copernicus Land Monitoring Service [60].

#### *2.4. Methods*

The proposed approach for monitoring post-fire vegetation state and estimating its dynamics included the following basic steps. First, spectral indices for the period between 2016 and 2021 were calculated to assess forest vegetation recovery dynamics. Second, statistical regression analyses using three of the spectral indices as variables were performed. The indices involved in the linear regression analyses were DI, MCARI2, and MSI. These indices are representative of key post-fire characteristics of the affected territories. The DI was used for the assessment of disturbance of forest ecosystems, burn severity, and vegetation damage. MCARI2 is representative for vegetation regrowth, and MSI shows stress in ecosystems caused by moisture deficiency. The third step considered differences in landscapes and the conditions for forest recovery they create, which is predetermined by the impact of slope exposures and their influence on the heat–moisture ratio. The assessment of the slope exposure factor was based on a differentiated evaluation of indices dynamics and their interpretation. The final step consisted of the validation of obtained results using statistical regression analyses involving the forest HRLs as independent reference data. Interactive three-dimensional panoramas from Google Street View were also used in the validation process. The interactive panoramas were used to generate photographs in different directions in order to demonstrate the state of various ecosystems. X and Y coordinates and altitude of the point locations of the photographs were extracted. The obtained point locations were georeferenced and digitized to be used in overlay analyses with the indices' rasters. The eight locations affected by the fire that were observed via this technology were linked to the obtained results as a means for field observation. For each location, indices values were extracted by taking into account the observed perspective and the distance of the objects from the point of capture.

Spectral Indices Selected for Assessment of Post-Fire Vegetation Recovery

The spectral indices presented in Table 1 were selected and calculated to assess the post-fire vegetation state.

The most well-known and used vegetation index for quantifying green vegetation in the near-infrared wavelength region and chlorophyll absorption in the red wavelength region is the NDVI. The NDVI strongly correlates with climate variations and their impact on plant growth. That makes this index especially suitable for estimations of climate-related vegetation changes. Moreover, changes in NDVI values correlate with the de Martonne aridity index [19].

The MCARI has been proposed to estimate variations in chlorophyll content and its concentration changes. Unlike MCARI, the newly designed MCARI2 is less sensitive to chlorophyll concentration variations but has a high linear relationship with near-infrared canopy reflectance and high linearity with the green LAI. The LAI is an important variable used to estimate the biophysical processes of different vegetation types and predict their growth and productivity [21,62]. Both the NDVI and MCARI2 range between −1 and +1. The highest values indicate dense and "healthy" vegetation, and the lowest values indicate dead plants or inanimate objects.


**Table 1.** Spectral indices calculated in this study.

The remotely sensed MSI is used for canopy stress analysis, productivity prediction, and biophysical modeling. It detects plant water stress for these plants only, which are able to tolerate low leaf water content through cellular adjustments. In the current study, the MSI is especially suitable for monitoring coniferous forests and assessing specific damages that cannot be detected using NIR/R vegetation indices [39]. Considering coniferous vegetation, the differences in the MSI between damaged and undamaged stands are not necessarily related to differences in the LAI. The MSI ranges from 0 to more than 3. Higher values indicate greater moisture stress.

Additionally, a quantitative analysis of forest degradation resulting from forest fires was performed by introducing the NDGI and the NDWNI. Both indices are based on satellite image orthogonalization. In the process of orthogonalization, three differentiable classes (soil brightness, greenness, and wetness axes) related to the main components of the Earth's surface (soil, vegetation, and water) are obtained. The NDGI uses the greenness component, which corresponds to vegetation's spectral reflectance characteristics (SRC), and the NDWNI uses the wetness component, which corresponds to water's SRC. These indices quantitatively estimate the slightly positive and negative values of change in the vegetation's green mass and moisture content for a given period [38,64,65]. Both indices range from −1 to +1. The positive NDGI values indicate plant growth and improvement of the vegetation state, and the negative values indicate deterioration of the vegetation state or deforestation. The positive NDWNI values indicate an increase in moisture content in ecosystems, and the negative values indicate a decrease in moisture content.

The DI has also been used to monitor disturbances by forest fires. The index values range widely, with positive values indicating disturbances. Higher values indicate more severe disturbance. The DI is modified by weighting each input component to maximize the difference between disturbed and undisturbed forest canopy. The weights reduce the effects of background variations while emphasizing the variations caused by disturbance [48]. The model for calculating the DI includes three steps: the first step is the decomposition of each of the three major Tasseled Cap components (brightness (BR), greenness (GR), and wetness (W)); the second step is to calculate the averages and standard deviations for each of the Tasseled Cap components; and the third step is to calculate the normalized values of the components. These steps are needed to normalize the radiometric changes. In the normalization, the following equations were used [46]:

$$\text{InBR} = \text{(BR} - \text{E\{BR\}} / \text{St.Dev (BR)} \tag{9}$$

$$\text{InGR} = \text{(GR} - \text{E} \text{\textasciicircum} \text{(\textasciicircum)} \text{(\text{St.Dev (GR)} \text{)}} \tag{10}$$

$$\text{InW} = (\text{W} - \text{E} \text{\textdegree\text{W}}) / \text{St.Dev \text{\textdegree}} \text{ (W)} \tag{11}$$

where E {BR}, E {GR}, and E {W} are the average values of the brightness, greenness, and wetness, respectively. St.Dev (BR), St.Dev (GR), and St.Dev (W) are the respective standard deviations of these Tasseled Cap components. Therefore, nBR, nGR, and nW are the normalized values of brightness, greenness, and wetness, respectively. The DI is computed according to the equation presented in Table 1.

#### **3. Results**

#### *3.1. Forest Vegetation Recovery for the Entire Study Area*

Areas with negative NDGI values predominated until 2019. The NDGI values were negative even a year before the fire event, which indicates a degraded state of forest vegetation. The territories with negative NDGI values had maximum territorial spread in the pre-fire year and a year after the fire (Figure 3, Table A2). The areas in the range from −1 to −0.8 were dominant between 2016 and 2017 in the year immediately following the fire (Figure 3, Table A2). With increasing temporal distance from the fire event, the maximum territorial spread shifted to the category with slightly positive values (from 0 to 0.2). The areas with positive NDGI values, or an increase in biomass, had minimal spread in the year immediately following the fire (Figure 3, Table A2). The positive NDGI values had the highest growth between 2019 and 2020 (Figure 3, Table A2).

Surprisingly, the areas with no disturbance or negative DI values prevailed during the entire studied period (Figure 3, Table A2). However, this maximum was highest in the year before the fire event. Moreover, during the same year, the highest territorial spread of the areas with great disturbance was observed. Areas from the same category (DI > 5) were also recorded in the following 2017 and 2018 years (Figure 3, Table A2). After 2019, with increasing temporal distance from the fire event, areas falling into this category were not observed. In the year before the fire event, a large portion of the areas had DI values between 0 and 2. In 2017, the areas were almost evenly distributed in the categories with DI values between 0 and 4. In the following four years, between 2018 and 2021, the maximum spread had areas with DI values between 0 and 3 (Figure 3, Table A2).

The first post-fire year was characterized by the most pronounced stress due to moisture deficiency. In the same year, the highest territorial spread of the category with MSI values above 1.5 was also observed (Figure 3, Table A2). The MSI category with values between 0.5 and 1.5 had the largest share of the area. The maximum spread had territories with MSI between 1.1 and 1.3. In the pre-fire year, 60.5% were in the category with MSI between 0.5 and 0.7. Between 2018 and 2021, the area was concentrated in the MSI categories between 0.5 and 1.3. The maximum spread was between 0.7 and 1.3 (Figure 3, Table A2).

Regarding the NDWNI, the maximum area was mainly in the category between −0.2 and 0.2, indicating weak dynamics in the moisture content. An exception was the oneyear period immediately after the fire when a significant increase in moisture content was observed (Figure 3, Table A2). In the years before the fire and between 2018 and 2021, a significantly smaller share of the area had positive NDWNI dynamics (Figure 3, Table A2).

During the year before the fire, 94% of the area had MCARI2 values above 0.6. Moreover, almost half of the area had MCARI2 values between 0.7 and 0.8. This category (0.7–0.8) also had a maximum territorial spread in the post-fire years, but the area falling into this category was significantly less (Figure 3, Table A2). In the post-fire year of 2017, the area was more evenly distributed between the individual categories. In the following years, the area was concentrated mainly in the MCARI2 category with values between 0.4 and 0.9. The share of the areas in this category gradually increased, reaching 96% in 2021 (Figure 3, Table A2).

**Figure 3.** Dynamics of the spectral indices calculated for the period of 2016–2021.

Greater dynamics were observed in the maximum territorial spread between the individual NDVI categories. In the pre-fire year, the highest proportion of territories had NDVI values between 0.6 and 0.7, whereas in the year following the fire event, maximum spread had territories with NDVI between 0.3 and 0.4. The maximum territorial spread shifted to the categories with higher NDVI values in the years between 2018 and 2021 (Figure 3, Table A2).

#### *3.2. Forest Vegetation Recovery in the Individual Slope Aspects*

The differential analysis of the dynamics of indices values on the individual slope exposures showed that in the first three years after the fire event (2017–2019), the southwestfacing slopes had faster vegetation recovery and less moisture stress. These slopes had the highest values of the NDVI and MCARI2 and the lowest values of the MSI (Table 2). The MCARI2 had a better ability to differentiate vegetation recovery through the individual

slopes. On east- and southeast-facing slopes, the NDVI showed equal totals over the three years, while the MCARI2 distinguished them. The difference in MCARI2 values between the individual slope aspects was higher than those of the NDVI and was clearer (Table 2). Regarding moisture stress in the first three years after the fire event, the MSI total values gradually increased in the following sequence: northeast-facing, southeast-facing, eastfacing, and south-facing slopes (Table 2). In the last two post-fire years, the northeast-facing slopes had the highest total values of NDVI and MCARI2 and the lowest moisture stress. On northeast-facing slopes, the highest totals of NDVI and MCARI2 and the lowest totals of MSI for the entire period were recorded (Table 2).


**Table 2.** Mean values of the spectral indices for the different slope aspects in the studied years.

In the pre-fire year, only the vegetation on southeast-facing slopes had positive DI values (Table 2). Moreover, territories most affected by the fire were situated on slopes with such exposure in the northern part of the area. The entire forest vegetation was destroyed in those territories (Figure 1). The fire most likely started there (Figure 3). The mean DI values were positive for all slopes in the following three post-fire years. The total DI values for the individual slope exposers for 2017, 2018, and 2019 decreased in the following sequence: south-facing, east-facing, northeast-facing, southeast-facing, and southwest-facing slopes (Table 2). Despite the higher DI values recorded for the east and northeast-facing slopes in the first three post-fire years, only these slopes had negative values in the last two years of observation (Table 2). On slopes with a southern component, the DI remained positive. In contrast to the trend observed in the NDVI, MCARI2, and MSI for the entire period, the lowest DI values and the optimal vegetation state were not recorded for the northeastern but for the southwestern slopes. However, as recorded by the listed indices, the vegetation had the worst condition on south-facing slopes (Table 2).

The other two indices based on TCT (NDWNI and NDGI) can measure small changes in moisture content and green mass. A slight positive increase in moisture content in the first post-fire year was recorded only for the southeast-facing slopes. In 2018 and 2019, the northeast-facing slopes had the highest positive NDWNI dynamics. The eastern slopes also showed positive dynamics in the moisture content in these years (Table 2).

The NDGI, or green mass, showed positive dynamics after the second post-fire year (2017–2018). In the period between 2017 and 2018, the vegetation on east- and northeastfacing slopes had the highest increase in NDGI values. These were the highest NDGI mean values in the entire period of observation and for all of the slope exposures (Table 2). Moreover, during the drought in 2020, the vegetation on the northeastern and eastern slopes again showed better condition compared to those on the slopes with different exposure (Table 2). The NDWNI dynamics were also positive for the northeast- and east-facing slopes in 2020. However, for the entire period of observation, the vegetation on southeast-facing slopes had the highest total increase in NDGI values (Table 2).

#### *3.3. Linear Regression Analyses for the Apectral Indices*

The statistical regression analyses were representative for the correlation and dependency between some of the key post-fire characteristics of the affected territories in the vegetation recovery process.

The correlation between the areas affected by disturbance increased with the development of ecosystem restoration processes and improving vegetation condition. The weaker correlation between 2017 and 2018 indicated a higher difference in the state of the forest vegetation in the first and second post-fire years (Figure 4c). This trend confirmed the rapid development of the post-fire succession process between the first and second years after the wildfire. The correlation was also lower between 2019 and 2020 when the severe drought in the summer of 2020 influenced the vegetation state and increased the intra-territorial differences. The increase in correlation for the periods between 2018 and 2019 and between 2020 and 2021, in turn, showed greater similarity in the state of the vegetation. In the second post-fire year (the period between 2018 and 2019), the higher correlation indicated delaying of the ecosystem restoration process, and in the 2020–2021 period, it indicated restoring balance in the ecosystems disturbed by the severe drought. The character of the territorial spread of vegetation throughout the individual categories divided by the DI values confirmed these observations (Figure 3, Table A2). In the first post-fire year, which was distinguished by weaker correlation, the vegetation was divided between eight categories. In 2021, as vegetation recovery progressed and the highest correlation was recorded, the number of these categories decreased to five. In 2021, 67.4% of the area was concentrated in two categories. 30% of this share consisted of the territories with DI between 1 and 2. That was the greatest share of the area falling into this category of disturbance for the entire period of observation. The area in this category gradually increased between 2016 and 2021, starting at 16.4% in 2016 and reaching up to 30% in 2021 (Figure 3, Table A2).

**Figure 4.** Linear regression analyses between DI and MCARI2 (**a**); DI and MSI (**b**); and DIs for two consecutive years (**c**) for the period between 2017 and 2021.

The correlation between MSI and DI values, in turn, showed a decreasing trend with the progression of the ecosystem restoration process (Figure 4b). In 2017 (the first post-fire year), the correlation between the disturbed ecosystems and moisture stress was highest for the entire period of observation. That indicated that for most of the forest ecosystems, higher MSI values were connected with higher DI values. In other words, the lack of

moisture was connected with a higher degree of disturbance of the ecosystems (Figure 4b). The progression of ecosystem restoration processes decreased this dependency.

The correlation and the dependency between MCARI2 and DI (i.e., between leaf area and intensification of photosynthesis and the degree of disturbance of the ecosystem) showed a decreasing trend with the progression of ecosystem restoration processes for the first three post-fire years. Figure 4a shows that the areas characterized by high MCARI2 values were distinguished by lack of disturbance of the ecosystem and vice versa; areas with lower MCARI2 values (e.g., 0.2) (i.e., areas with underdeveloped vegetation and small leaf area) were distinguished by a higher degree of disturbance of the ecosystems. The progression of ecosystem restoration processes led to an increase in leaf mass and intensification of photosynthesis. The lowest MCARI2 values reached 0.4 in 2019. In these areas, the DI had decreased to below 4. The decline in the dependency between the vegetation state and the degree of disturbance of ecosystems confirmed the progression of ecosystem restoration processes. The decreasing trend in the correlation was interrupted by drought stress in 2020. In 2021, the correlations between both indices started to decline again. However, the forest recovery trend still could be observed. Starting from 2020, vegetation with DI values above 4 was not observed (Figure 4a).

#### *3.4. Validation through HRLs*

Statistical regression analyses involving the forest HRLs acquired in 2012, 2015, and 2018 as independent reference data were performed to validate the obtained results.

Generally, between 2012 and 2018, an increase in tree cover density of both broadleaved forests and coniferous forests was observed. The largest share of the area had forests with tree cover density between 40% and 80% (Table 3). In the post-fire year of 2018, the non-forested areas, as expected, had the highest disturbance (Table 4). They had positive DI values. The DI for the forested areas in the same year was negative. The non-forested area recorded the highest moisture stress (MSI > 1) (Table 4) but also the highest increase in NDGI values. That was related to the development of vegetation succession processes in the deforested territories. The vegetation indices NDVI and MCARI2 recorded their highest values in broad-leaved forests. However, coniferous forests had the highest positive dynamics in moisture content (NDWNI) (Table 4).

**TCD (%) Broad-Leaved Forests Coniferous Forests 2012 (%) 2015 (%) 2018 (%) 2012 (%) 2015 (%) 2018 (%)** 0–20 6.66 0.8 0.72 0.63 non non 20–40 14.53 9.41 25.51 3.83 0.75 9.39 40–60 37.61 72.47 39.02 19.15 43.85 54.26 60–80 38.97 15.92 28.47 73.96 55.35 34.96 >80 2.22 1.4 6.29 2.43 0.06 1.39 Mean 43.86 52.21 56.19 47.37 55.67 60.28

**Table 3.** Tree-cover density (TCD) classes (%) and area of their distribution (%) within broad-leaved and coniferous forests in 2012, 2015, and 2018.

The DI and MSI were the indices that showed the highest correlation with forest density (Table 5). The correlation between the DI and forest density was the highest in the coniferous forests. In the broad-leaved forests, the correlation between these two variables was slightly lower. Regarding the moisture content, both deciduous and coniferous forests showed a similar correlation. Generally, the coniferous forests were distinguished by higher dependence on the spectral indices. The difference in correlation between the forest density and NDVI was significant when comparing deciduous and coniferous forests. In deciduous forests, R was barely 0.162, whereas in coniferous forests, it was 0.529 (Table 5). The observations for the NDGI were similar. The difference between both forest types

was greater than four times. The correlation between forest density and MCARI2 was also slightly higher in the coniferous compared with broad-leaved forests (Table 5).


**Table 4.** Mean values of the spectral indices in broad-leaved forests, coniferous forests, and in non-forested areas in 2018.

**Table 5.** Linear regression between forest density (independent) and spectral indices (dependent) for each of the forest types.


#### *3.5. Validation through Field Observation*

Interactive three-dimensional panoramas of eight-point locations from Google Street View were also used in the validation process. The eight locations were linked to the obtained results as a means for field observation. The indices values clearly showed the trends in the manifestation of the indices in relation to the observed ecosystems' components (i.e., soil and vegetation). The lowest indices values were observed in the post-fire karst territories in point one and point eight, and the highest values were observed in the forest territories unaffected by the fire (Figure A1, Table A3).

#### **4. Discussion**

The deteriorated vegetation state and the landscape-ecological conditions (karst terrain, dry and hot summers, and warm slope exposures) caused the fire. Dry vegetation and steep terrain made the fire more intense, devastating, and hard to control. The consequences in the area of occurrence were almost complete destruction of the forest vegetation, litter cover, and soil organic layer.

The soil type and slope exposures are the main landscape-forming factors that determine differences in the processes of vegetation recovery in the study area. Unfavorable characteristics of the soils in this area, such as shallow profile, low organic matter content, acidic soil reaction, and low exchange capacity, worsened after the fire and further inhibited the recovery of vegetation, including forest. These processes manifested more significantly on the southern slopes. The removal of vegetation and litter cover as a result of a fire reduces rainfall interception, which enhances runoff and erosion rates [66]. Moreover, the burn severity of the soil surface on the south-facing slopes decreases soil carbon content and changes soil acidity [31].

The neighboring territories were also affected to a great extent. Forest management practices include a sanitary logging of burnt forest stands after a fire. For this reason, in 2018 (two years after the fire), actions to remove the burnt forest vegetation were taken. As a result, a large part of the territory was completely cut down. The recovery processes of the forest vegetation were interrupted. The manifestation of TCT-based indices (NDGI and NDWNI), which represent the dynamics in greenness and wetness after 2017–2018, was also an indication for this interruption of the recovery process (Figure 3). It was evident from a sharp decline of the areas with values that ranged between −1 and −0.8, which is typical for forested territories, and from a significant increase of the areas in the NDGI category of 0.8–1 (Figure 3, Table A2), which indicates the rapid development of grass vegetation.

The negative values of the TCT-based DI were also indicative of the landscapes' dynamic conditions in the territories affected by the fire. This index showed stable presence of areas "undisturbed" by the fire. Their share was most significant (52% of the studied area) before the fire (10 July 2016). In the years after the fire, their share remained relatively stable at about 1/3 of the territory (Figure 3, Table A2). The analysis showed that these were forest areas that remained undamaged by the fire as well as local patches of meadows in the southern zone unaffected by the fire. This category also included areas with exposed bedrock formed due to the nature of lithological type (karst) in the northern part of the study area [54]. In these areas, the fire was not a destabilizing factor for the ecosystems' condition (Figure 1).

The maximum spread of the areas in the NDWNI category (−0.2–+0.2), which had a percentage distribution similar to that observed for the negative values of DI (approximately 30%) of the studied area (Figure 3, Table A2), also confirmed the observations made in the analysis of the DI. The territories unaffected by the fire, were distinguished by the lowest dynamics in moisture content. The first post-fire year was distinguished by the highest NDWNI dynamics. It was related to rapid post-fire succession processes in a large part of the study area. However, the NDGI values, which were representative of greenness dynamics, did not indicate high dynamics in the same year, but they did indicate high dynamics in the following year after sanitary logging (Figure 3, Table A2). The increase in moisture content in the first year after the fire was not induced primarily by the better state of vegetation but was mainly a result of the meteorological condition in the summer of 2017 when the precipitation was 226% more than what it typical for the region (Figure 2b). It should be noted that in the first post-fire years, vegetation primarily comprises herbaceous species. The state of this type of vegetation is closely dependent on environmental conditions. Grasslands are less resistant to anomalies related to temperature and humidity [67,68]. Moreover, they are strongly dependent on their fluctuations [69].

The indices that were not based on orthogonalization showed the general dynamics in the vegetation state. Generally, the areas with comparatively high values of NDVI and MCARI2 (taking into consideration the character of vegetation), which corresponds with lower MSI values, had the greatest territorial spread in the pre-fire year (Figure 3, Table A2). In the first post-fire year, the areas were distributed among a large number of categories, i.e., the vegetation in the fire-affected territories was characterized by great diversity in terms of its state. There were both unaffected areas and areas with varying degrees of post-fire damage. With an increasing temporal distance from the fire event, an increasing share of the territories were concentrated in fewer and fewer categories, with a slight shift towards the categories representing a better state of vegetation. In the different post-fire years, there was a slight growth or contraction of the individual categories within the general trend of vegetation recovery (Figure 3, Table A2). However, this dynamic was the result of climatic elements in the relevant year of observation (Figure 2a,b).

Some trends stand out regarding the influence of slope exposures on the vegetation recovery process. Vegetation recovery was faster on warmer slopes in the first post-fire years (2017–2019) (Table 2). However, it should be noted that during those years, the lack of moisture was not a limiting factor for the recovery processes (Figure 2b). In all three years, the summer precipitation exceeded the norm (Figure 2b). These results are consistent with those of Wilson at al. [70], who found that vegetation had a higher recovery rate when the temperature and precipitation were higher. A rapid, "catch-up" development on the northeastern slopes was observed (Table 2) in the last two years of observation (2020 and 2021). These were dry years (the precipitation sum for the three summer months was below the norm), and the drought delayed vegetation recovery on the warmer slopes (Figure 2b). The recovery processes on the northeastern slopes for the last two years of observation were so rapid that they influenced the values indicating the total recovery for the entire period (Table 2). This rapid vegetation recovery was also recorded by the DI. The tendency of positive values (i.e., values indicating disturbance in ecosystems) was interrupted only for the slopes with northeast and east exposure (Table 2). The DI showed a significant relationship with post-fire vegetation recovery processes and that its dynamics, when assessing such a process, were correlated with various climatic and topographic factors [42]. Chen et al. [42] confirmed that the role of slope exposures in the dynamics of post-fire vegetation recovery was essential. Our results are consistent with those of Chen et al. [42]. They found that vegetation recovery on sunny sides was greater than on shady sides and that it closely related to elevation and its influence on the heat–moisture ratio [42].

The NDWNI and NDGI showed the highest positive dynamics for the north- and northeast-facing slopes. These slopes are in the western part of the study area, where due to the lower degree of slope inclinations, the forest vegetation was largely preserved from the fire (Figure 1). After the post-fire sanitary logging in 2018, small patches of these slopes were deforested. These patches remained between separate groups of trees (Figure 1). Regarding the landscape-forming factors (soil type [54], slope aspects, and inclinations), these territories had more favorable conditions for vegetation recovery.

The results of the statistical and validation analyses indicated the reliability of the methodology used to monitor and assess post-fire vegetation recovery processes. The linear regression analyses showed stronger correlations between the objects of the Earth's surface that were under favorable conditions for vegetation recovery for two consecutive years (higher R values between DI) and vice versa (weaker correlations when the territory was under post-fire or drought stress in one of two consecutive years (lower R value between DI)) (Figure 4c). The dependency between the severity of disturbance and moisture content decreased with increasing time from the fire event and with the development of the vegetation recovery process (Figure 4b). The correlation between the leaf area (indirectly represented by MCARI2) and the disturbance of the ecosystems decreased with the development of vegetation recovery. This trend of decreasing correlation was interrupted by the drought in 2020 that impacted the general state of the vegetation and resulted in withering and shrinking of tree leaves (Figure 4a).

The validation through high-resolution forest layers showed higher dependency between the forest density and the indices values for the coniferous forests (Table 5). This tendency is induced by the fact that coniferous vegetation, in contrast to deciduous vegetation, stays green during the entire year. Satellite images acquired in August or early September for each of the years were used. At the end of summer, deciduous forests undergo senescence. The process leads to reduced greenness and moisture in the tree leaves, disrupting the process of photosynthesis. In coniferous vegetation, these processes are significantly less noticeable. For that reason, the correlation between the forest density and the indices values representing the general state and functioning of vegetation were higher in the coniferous forests compared with the deciduous forests.

#### **5. Conclusions**

This study traced post-fire vegetation recovery dynamics using two groups of spectral vegetation indices and taking into consideration the local landscape factors in an area affected by a fire in southeast Bulgaria. Consistent with other studies [42,70], it was confirmed that vegetation recovery was dependent on climatic factors [70] and topography features [42]. Regarding the effectiveness of spectral vegetation indices for monitoring the post-fire vegetation state, it can be summarized that the statistical analysis and validation procedures confirmed their reliability for the assessment of restoration processes of vegetation after a fire. The NDVI, MCARI2, and MSI indicated general trends in post-fire vegetation dynamics, and the TCT-based indices (DI, NDGI, and NDWNI) were found to be

suitable for more precise analyses of intra-territorial differences. The DI was advantageous for the differentiation of post-fire severity in ecosystems. The obtained results clearly showed the intra-territorial heterogeneity of post-fire vegetation recovery and the influence of local environmental factors on the dynamics of the process. The study demonstrated the need of multi-factor analysis in post-fire monitoring and could serve as a basis for further post-fire-related studies. Estimations of the impact of soil erosion would be particularly valuable. In such a study, the changes in the soil characteristics and in-depth analysis of slope steepness must be taken into account.

**Author Contributions:** Conceptualization, D.A., E.V. and L.F.; methodology, D.A., E.V. and L.F.; software, D.A.; validation, D.A.; formal analysis, D.A.; investigation, D.A.; resources, D.A. and L.F.; data curation, D.A.; writing—original draft preparation, D.A.; writing—review and editing, D.A.; visualization, D.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Spectral (in microns) and spatial (in meters) resolution of Sentinel-2 sensor.


#### **Appendix B**

**Table A2.** Territorial spread (in %) of the individual categories divided for each of the spectral indices within the study area. This table represents the values behind the output rasters from Figure 3.


**Table A2.** *Cont*.


#### **Appendix C**

**Figure A1.** Point locations of field observations.

**Table A3.** Spectral indices values for each location. Interactive three-dimensional panoramas from Google Street view were used to generate photographs.



#### **Table A3.** *Cont*.

#### **References**

