Phenology Patterns and Postfire Vegetation Regeneration in the Chiquitania Region of Bolivia Using Sentinel-2

The natural regeneration of ecosystems impacted by fires is a high priority in Bolivia, and represents one of the country’s greatest environmental challenges. With the abundance of spatial data and access to improved technologies, it is critical to provide an effective method of analysis to evaluate changes in land use in the face of the global need to understand the dynamics of vegetation in regeneration processes. In this context, we evaluated the dynamics of natural regeneration through phenological patterns by measuring the maximal and minimal spectral thresholds at four fire-impacted sites in Chiquitania in 2019 and 2020, and compared them with unburned areas using harmonic fitted values of the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR). We used two-way ANOVA test to evaluate the significant differences in the values of the profiles of NDVI and NBR indices. We quantified severity at the four study sites using the dNBR obtained from the difference between preand postfire NBR. Additionally, we selected 66 sampling sites to apply the Composite Burn Index (CBI) methodology. Our results indicate that NBR is the most reliable index for interannual comparisons and determining changes in the phenological pattern, which allow for the detection of postfire regeneration. Fire severity levels based on dNBR and CBI indices are reliable methodologies that allow for determining the severity and dynamics of changes in postfire regeneration levels in forested and nonforested areas.


Introduction
Bolivia is among the countries with the highest wildfire activity recorded in South America in recent years [1,2]. In the last 20 years (2001-2020), wildfires in Bolivia directly affected an area of 23.6 million hectares, with the most extensive burnt areas being reported in 2004, 2005, 2010 and 2019 [3]. Annually, savannas are the most affected ecosystems by wildfires, followed by shrublands and forests [3]. Fire is mainly caused by human actions, principally due to pastures renewing the practice of slash and burn, known as chaqueo [4].
The Department of Santa Cruz, located in the eastern region of Bolivia, is one of the most affected by the wildfires, which mainly occur in the lowlands. In 2019, the worst occurrences of forest fires were recorded in Santa Cruz. As a result, 3.7 million hectares were affected, mainly in the Chiquitania region, of which 60% were forested areas [5]. These fires caused a series of ecological (e.g., impacts on soils, habitat losses, changes in species composition), economic (e.g., timber and wood products lost, agricultural production lost) and social (e.g., injured, air quality and smoke impacts) problems on the region, and directly affected the livelihoods of around 370 local communities [6]. In 2020, Santa Cruz recorded fires with 2.2 million hectares burned, of which 58% were forests [7]. These events reopened the debate about forest fires in the country, proposing different actions to recover Guasu Protected Area of Conservation and Ecological Importance (Ñembi Guasu, Figure  2). The criteria for selecting these sites were based on the presence of burned areas between 2019 and 2020, favorable conditions for entry permits by local populations, levels of legal protection, areas with productive agricultural or livestock activities, ease of access to sites impacted by fires, and ecosystem heterogeneity.

Data Sources and Methods
For the different analyses, we used Sentinel-2 satellite data that provide images from 10 to 60 m spatial resolution with 13 spectral bands and five-day revisit frequency [26]. We generated fire scar maps for 2019 and 2020 from Sentinel-2 Level 1-C reflectance satellite imagery (TOA) with cloud cover less than 20%, which are publicly available on the USGS EarthExplorer portal (https://earthexplorer.usgs.gov, accessed on 10 March 2021). We performed different analyses in Google Earth Engine (GEE) with Sentinel-2 images to obtain the spectral indices of Normalized Difference Vegetation Index (NDVI [27]) and Normalized Burn Ratio (NBR [28]), using QA60 band to mask clouds. The NDVI shows values between -1 and 1, and is determined using the near-infrared (NIR) and red (RED) bands (Table 2), where values close to -1 indicate little photosynthetic activity and thereby little growth of or reduction in vegetation, while values close to 1 reflect the opposite. The NBR is obtained from the NIR and shortwave infrared (SWIR) bands ( Table 2) that provides a better distinction between burned and unburned areas, as well as an optimal signal to obtain information on the variation of fire severity [28]. We used the land cover map (forest, nonforest, anthropogenic) elaborated by Maillard et al. [5] to intersect them with the 2019 and 2020 fire scar mapping, and thus determine the differences between burned and unburned natural areas. The obtained confidence level from the coverage classification was 85% [5]. Table 2. Spectral indices tested in the study, and their formulation and related spectral bands (band numbering refers to the Sentinel-2 sensor system).

Forest Fire Scars
We identified burned areas for 2019 and 2020 on the basis of the unsupervised classification of 56 Sentinel-2 scenes (Figure 1) for the July-December period using the Interactive Self-Organizing Data Analysis (ISODATA) algorithm [29] in ERDAS Image software to separate pixels into groups with similar spectral signatures. Previous studies that

Data Sources and Methods
For the different analyses, we used Sentinel-2 satellite data that provide images from 10 to 60 m spatial resolution with 13 spectral bands and five-day revisit frequency [26]. We generated fire scar maps for 2019 and 2020 from Sentinel-2 Level 1-C reflectance satellite imagery (TOA) with cloud cover less than 20%, which are publicly available on the USGS EarthExplorer portal (https://earthexplorer.usgs.gov, accessed on 10 March 2021. We performed different analyses in Google Earth Engine (GEE) with Sentinel-2 images to obtain the spectral indices of Normalized Difference Vegetation Index (NDVI [27]) and Normalized Burn Ratio (NBR [28]), using QA60 band to mask clouds. The NDVI shows values between -1 and 1, and is determined using the near-infrared (NIR) and red (RED) bands (Table 2), where values close to -1 indicate little photosynthetic activity and thereby little growth of or reduction in vegetation, while values close to 1 reflect the opposite. The NBR is obtained from the NIR and shortwave infrared (SWIR) bands ( Table 2) that provides a better distinction between burned and unburned areas, as well as an optimal signal to obtain information on the variation of fire severity [28]. We used the land cover map (forest, nonforest, anthropogenic) elaborated by Maillard et al. [5] to intersect them with the 2019 and 2020 fire scar mapping, and thus determine the differences between burned and unburned natural areas. The obtained confidence level from the coverage classification was 85% [5]. Table 2. Spectral indices tested in the study, and their formulation and related spectral bands (band numbering refers to the Sentinel-2 sensor system).

Forest Fire Scars
We identified burned areas for 2019 and 2020 on the basis of the unsupervised classification of 56 Sentinel-2 scenes (Figure 1) for the July-December period using the Interactive Self-Organizing Data Analysis (ISODATA) algorithm [29] in ERDAS Image software to separate pixels into groups with similar spectral signatures. Previous studies that identified forest fire scars in the Chiquitania region of Bolivia obtained accurate classification results using ISODATA [5,7]. Classification was performed with 100 maximal number of iterations, a convergence threshold of 100% (percentage of pixels whose class values should not change between iterations), and skip factors of 1. This resulted in a thematic classification of 25 spectral classes. A combination of image-specific RGB bands was also used for a visual review of the burned area detection and classification procedure. Obtained results were vectorized in ESRI ArcMap. To evaluate the level of uncertainty of the resulting classification, 250 field verification points were obtained between 2019 and 2021. In total, 66 sampling plots were obtained through the Composite Burn Index (CBI), 150 sampling points with ArcGisSurvey123 application, and 34 points with high-resolution images taken with UAV. This approach allowed for obtaining a confusion matrix [30] and three types of accuracy estimates: overall accuracy, user accuracy (commission error), and producer accuracy (omission error), including their 95% confidence intervals [30]. The overall accuracy obtained from the fire scars was 95.48% (Table 3).

Phenological Patterns of Burned and Unburned Areas
We evaluated the phenological patterns of vegetation in burned and unburned areas, identifying the thresholds of the NDVI and NBR indices. Both indices highlight changes in vegetation reflectivity when affected by fire. In GEE, we developed a script to perform probability sampling with random point spread, stratified by the cover type and proportion of burned and unburned areas for each study site. The sample number was calculated as proposed by Olofsson et al. [30], with a 90% confidence level for forests and 95% for nonforests. For the four study areas, a total of 1418 points were obtained, 393 for burned and 1025 for unburned areas. In GEE, we employed harmonic regression, a mathematical technique used to disaggregate a complex static signal into a series of individual sine or cosine waves, each characterized by a specific phase and amplitude [31]. We used the harmonic regression algorithms developed by Clinton [32] for Sentinel-2 collections to estimate the medians and harmonic fitted NDVI and NBR values for the selected land cover categories ( Figure 3). The amplitude of harmonic regression is indicative of the magnitude of changes in spectral bands seasonally, which could be related to lapsing. To classify the phenological patterns of the mainland cover categories at different phases, we conducted time series analysis of the NDVI and NBR indices. We then identified the phenological patterns (monthly, seasonal, and annual) of burned and unburned areas, with the definition of average annual maximal and minimal thresholds for the period of 2016-2020, excluding outliers.

Fire severity Levels Based on dNBR
Fire severity is an important measure in determining postfire vegetation resilience and succession [33]. We quantified severity at the four study sites using the differential NBR or Delta NBR index (dNBR or ∆NBR, delta normalized burn ratio) obtained from the difference between pre-and postfire NBR (Table 2). Higher dNBR values indicate more severe damage, while areas with negative dNBR values may indicate new vegetation growth after a fire. On the basis of Key and Benson [28], we classified dNBR results into seven severity levels: enhanced regrowth/low (<−0.25), enhanced regrowth/high (−0.25 to −0.1), unburned (−0.1 to 0.1), low severity (0.1 to 0.27), moderate-low severity (0.27 to 0.44), moderate-high severity (0.44 to 0.66), and high severity (>0.66). Although the dNBR category thresholds proposed by Key and Benson [28] were developed for boreal ecosystems, we used these classes because they come from a globally standard widely used methodology to assess fire severity.

Fire Severity Based on Composite Burn Index (CBI)
The Composite Burn Index [28] (CBI) is based on a visual assessment of the quantity of consumed fuel, the degree of soil charring, and the degree of tree mortality [34], and it was applied to and adjusted in a wide variety of ecosystems [35]. Between 2019 and 2021, we surveyed burn severity in the four study sites (Table 4) using the CBI methodology. We evaluated fire affectation on the surface soil substrate; grasses, shrubs, and trees less than 1 m tall; shrubs and trees 1-5 m tall; intermediate trees 10-25 cm trunk diameter and 8-20 m tall; and large trees receiving direct sunlight [28]. We hierarchically aggregated attribute scores for the five strata from understory to canopy and then averaged the obtained values. We selected 66 sampling sites within the four study areas to apply the CBI (Table 4) aiming for heterogeneity in vegetation types within the identified burn scars. At each site, we installed 400 m 2 sampling plots at different levels of burn severity. In each plot, we described the site characteristics accompanied by four supporting photographic records, georeferenced to the four cardinal directions. The thresholds of CBI severity levels employed in this study were: (a) unburned = 0-0.25, (b) low = 0.25-1.25, (c) moderate = 1.25-2.25, and (d) high = 2.25-3 [36].

Fire Severity Levels Based on dNBR
Fire severity is an important measure in determining postfire vegetation resilience and succession [33]. We quantified severity at the four study sites using the differential NBR or Delta NBR index (dNBR or ∆NBR, delta normalized burn ratio) obtained from the difference between pre-and postfire NBR ( Table 2). Higher dNBR values indicate more severe damage, while areas with negative dNBR values may indicate new vegetation growth after a fire. On the basis of Key and Benson [28], we classified dNBR results into seven severity levels: enhanced regrowth/low (<−0.25), enhanced regrowth/high (−0.25 to −0.1), unburned (−0.1 to 0.1), low severity (0.1 to 0.27), moderate-low severity (0.27 to 0.44), moderate-high severity (0.44 to 0.66), and high severity (>0.66). Although the dNBR category thresholds proposed by Key and Benson [28] were developed for boreal ecosystems, we used these classes because they come from a globally standard widely used methodology to assess fire severity.

Fire Severity Based on Composite Burn Index (CBI)
The Composite Burn Index [28] (CBI) is based on a visual assessment of the quantity of consumed fuel, the degree of soil charring, and the degree of tree mortality [34], and it was applied to and adjusted in a wide variety of ecosystems [35]. Between 2019 and 2021, we surveyed burn severity in the four study sites (Table 4) using the CBI methodology. We evaluated fire affectation on the surface soil substrate; grasses, shrubs, and trees less than 1 m tall; shrubs and trees 1-5 m tall; intermediate trees 10-25 cm trunk diameter and 8-20 m tall; and large trees receiving direct sunlight [28]. We hierarchically aggregated attribute scores for the five strata from understory to canopy and then averaged the obtained values. We selected 66 sampling sites within the four study areas to apply the CBI (Table 4) aiming for heterogeneity in vegetation types within the identified burn scars. At each site, we installed 400 m 2 sampling plots at different levels of burn severity. In each plot, we described the site characteristics accompanied by four supporting photographic records, georeferenced to the four cardinal directions. The thresholds of CBI severity levels employed in this study were: (a) unburned = 0-0.25, (b) low = 0.25-1.25, (c) moderate = 1.25-2.25, and (d) high = 2.25-3 [36].

Statistical Analysis
We used a two-way ANOVA test to evaluate significant differences in the values of the profiles fitted with the harmonic model of the Sentinel-2 NDVI and NBR indices for the five-year analytical period (2016-2020) in relation to burned versus unburned areas for the four study sites in the Department of Santa Cruz. The results of two-way ANOVA are shown in Table 5. To determine the relationship between the CBI and dNBR indices for 2019, Spearman's correlation analysis was performed (p < 0.05). We used RStudio [37] for all statistical analyses.

Phenological Patterns of Burned and Unburned Areas
Analysis of the phenological patterns based on the fitted values of the harmonic model for Sentinel-2 NDVI and NBR showed marked seasonal (rainy season and dry season) and interannual variability (Figure 4) in the forested and nonforested areas of three of the four sites. The highest photosynthetic activity was recorded in the months of January-May and November-December, while June-October showed low chlorophyll absorption ( Figure S1). In the two-way ANOVA test, significant differences were evident, mainly for NBR, when comparing years, and burned and unburned areas (Table 5). These results indicate that NBR is the most reliable index for interannual comparisons and determining changes in the phenological pattern ( Figure 4). The main changes in NBR were between the minimal threshold levels of the prefire period (2016-2018) and the postfire values (2019-2020).
At AltaVista, interannual comparisons indicate a change in the minimal threshold of NBR in the forest that was impacted by fire, where the prefire stage recorded an average value of 0.31, which decreased to 0.16 in 2019 and reached 0.32 in 2020 (Figure 4), demonstrating a regeneration process one year after the fire. At Copaibo, the average minimal NBR threshold of the burned forest in the prefire stage was 0.61, which increased to 0.54 in 2019 and 0.34 in 2020. According to our analysis, no regeneration processes were detected in Copaibo due to the continuous fires in 2020 (Figure 4). At Laguna Marfil, interannual comparisons show some changes in the phenological patterns of the forest, with an increase in the minimal threshold from 0.34 in 2016-2018 to 0.22 in 2020 (Figure 4). It is evident that there is a regeneration process in this forest. However, in the nonforested areas of Marfil, the average values during the pre-and postfire periods were 0.14 and -0.15, respectively, indicating no regeneration (Figure 4). This was due to registered fires in two consecutive years (2019-2020). In Ñembi Guasu, the burned forest shows changes in the minimal NBR threshold, from a value of 0.31 in the prefire period to −0.43 (2019) and 0.17 (2020) in the postfire period (Figure 4). In the nonforested areas (Abayoy) of Ñembi Guasu that had been impacted by fire, the average value of the minimal NBR threshold was 0.25 (prefire), which increased after the fire event to −0.36 in 2019 and 0.17 in 2020. A regeneration process is evident in both types of Ñembi Guasu ecosystems.

Fire Severity Levels Based on dNBR
Values obtained in the dNBR for fire scars from the four study areas indicate a natural regeneration scenario between 2019 and 2020, even though some severity levels were not reduced because some sites burned in 2019 and 2020 ( Figures 5-8, Table S1). At Alta Vista, the fire categories with the highest proportion in 2019 were low severity (46%) and unburned (27%); however, by 2020, there was an increase in dNBR for the lower-severity categories, such as enhanced regrowth/high with 64% and enhanced regrowth/low with 34% ( Figure 5, Table S1). At Copaibo, fire severity levels in 2019 were mainly concentrated in the categories of moderate-low severity with 37%, low severity with 27%, and moderatehigh severity with 24%. Although a reduction in severity levels was expected, fires in 2020 caused these two categories to continue to present the highest percentages with 21%, 24% and 22%, respectively ( Figure 6, Table S1). At Laguna Marfil, the dNBR levels for 2019 and 2020 indicate a reduction in severity levels in moderate-low severity from 51% to 11%, and low severity from 25% to 17%; in 2020, the unburned and enhanced regrowth/high classes increased to 24% and 26%, respectively (Figure 7, Table S1). At Ñembi Guasu, the severity levels of fires in 2019 based on dNBR were mainly concentrated in the moderatehigh-severity category with 32%, and in the categories of moderate-low severity with 20% and moderate-high severity with 18%. By 2020, the low-severity (40%) and unburned (25%) categories increased considerably, demonstrating the recovery of vegetation in this protected area ( Figure 8, Table S1).

Fire Severity Based on CBI
Of the 66 CBI plots that we installed in the four study areas between 2019 and 2021, the 16 that we remeasured showed a natural regeneration process (Table S2). At Alta Vista, the three CBI plots in the Chiquitano Forest demonstrated moderate severity levels in 2019 (CBI = 1.4-2), and subsequently resulted in reduced severities of low and unburned categories in 2020 (CBI = 0.1-0.8) (Table S2). At Copaibo, 3 of the 14 severity plots that we measured in 2021 were not affected by the 2019 or 2020 fires (Table S2), so the CBI values were very low (<0.4); the remaining 11 showed severity levels in the unburned-to-moderate categories (CBI = 0.4-2.1).
At Ñembi Guasu, the 20 established CBI plots were located in areas impacted by fires in 2019 (Table S2). For the scrubland formation known as Abayoy, the field evaluations that we had conducted months after the 2019 fire series showed high severity levels (CBI = 2.3-2.7). In the Chiquitano Forest transition to Chaco of Ñembi Guasu (Table S2), fire severities also showed high levels (CBI = 2.4-2.7). Of the total sampling in the four study areas, both ecosystems presented the highest levels of severity.  (Figure 4). In the nonforested areas (Abayoy) of Ñembi Guasu that had been impacted by fire, the average value of the minimal NBR threshold was 0.25 (prefire), which increased after the fire event to −0.36 in 2019 and 0.17 in 2020. A regeneration process is evident in both types of Ñembi Guasu ecosystems.   growth/high classes increased to 24% and 26%, respectively (Figure 7, Table S1). At Ñembi Guasu, the severity levels of fires in 2019 based on dNBR were mainly concentrated in the moderate-high-severity category with 32%, and in the categories of moderate-low severity with 20% and moderate-high severity with 18%. By 2020, the low-severity (40%) and unburned (25%) categories increased considerably, demonstrating the recovery of vegetation in this protected area (Figure 8, Table S1).

Fire Severity Based on CBI
Of the 66 CBI plots that we installed in the four study areas between 2019 and 2021, the 16 that we remeasured showed a natural regeneration process (Table S2). At Alta Vista, the three CBI plots in the Chiquitano Forest demonstrated moderate severity levels  Relationship analyses between the dNBR and CBI indices for burned forested and nonforested areas, show a significant and high correlation for 2019 (R = 0.82, p < 0.001) (Figure 9). This indicates that the results obtained with the dNBR can be considered reliable. severities also showed high levels (CBI = 2.4-2.7). Of the total sampling in the four study areas, both ecosystems presented the highest levels of severity. Relationship analyses between the dNBR and CBI indices for burned forested and nonforested areas, show a significant and high correlation for 2019 (R = 0.82, p < 0.001) (Figure 9). This indicates that the results obtained with the dNBR can be considered reliable.

Discussion
Wildfires reduce the reflectance of visible-to-NIR wavelengths due to the loss of photosynthetic vegetation, and increase the reflectance of SWIR wavelengths due to increased char and ash [38,39]. The postfire vegetation regeneration processes had the opposite effect on spectral reflectance. In this sense, we calculated and evaluated the NDVI and NBR indices, and compared the results for burned and unburned areas, from forested and nonforested ecosystems, allowing for us to determine the index that offers the best capabilities to reflect postfire regeneration in the ecosystems of the Chiquitania region. The forested cover of the four study areas had a much higher NDVI and NBR than those of nonforested areas, which decrease during the dry season. However, annual differences were recorded for the period 2016-2020, between the averages of the maximal and minimal thresholds in unburned areas. These differences are related to a response of the vegetation to annual seasonal characteristics, such as temperature, precipitation, and humidity. In months with low photosynthetic activity, the vegetation enters a dormancy process to save energy during the dry season, as had been recorded in other deciduous forests in Bolivia [21]. It is important to consider these results to avoid errors in the interpretation of spectral thresholds in the comparison of burned and unburned areas. In addition, we must consider meteorological drought events that were recorded in recent years on a continental scale [40], which were also perceived and recorded in Bolivia [5]. Our results show that the NBR showed a better spectral response, as it allowed for us to identify significant differences between years, and between burned and unburned areas because the red/NIR-based vegetation indices are mainly related to the density of green leaves in a given area, and are thereby a good indicator of plant cover and vitality [41]. The SWIR-based vegetation indices are sensitive to the detection of dead or nonphotosynthetic vegetation in a postfire

Discussion
Wildfires reduce the reflectance of visible-to-NIR wavelengths due to the loss of photosynthetic vegetation, and increase the reflectance of SWIR wavelengths due to increased char and ash [38,39]. The postfire vegetation regeneration processes had the opposite effect on spectral reflectance. In this sense, we calculated and evaluated the NDVI and NBR indices, and compared the results for burned and unburned areas, from forested and nonforested ecosystems, allowing for us to determine the index that offers the best capabilities to reflect postfire regeneration in the ecosystems of the Chiquitania region. The forested cover of the four study areas had a much higher NDVI and NBR than those of nonforested areas, which decrease during the dry season. However, annual differences were recorded for the period 2016-2020, between the averages of the maximal and minimal thresholds in unburned areas. These differences are related to a response of the vegetation to annual seasonal characteristics, such as temperature, precipitation, and humidity. In months with low photosynthetic activity, the vegetation enters a dormancy process to save energy during the dry season, as had been recorded in other deciduous forests in Bolivia [21]. It is important to consider these results to avoid errors in the interpretation of spectral thresholds in the comparison of burned and unburned areas. In addition, we must consider meteorological drought events that were recorded in recent years on a continental scale [40], which were also perceived and recorded in Bolivia [5]. Our results show that the NBR showed a better spectral response, as it allowed for us to identify significant differences between years, and between burned and unburned areas because the red/NIR-based vegetation indices are mainly related to the density of green leaves in a given area, and are thereby a good indicator of plant cover and vitality [41]. The SWIR-based vegetation indices are sensitive to the detection of dead or nonphotosynthetic vegetation in a postfire environment [42]. In other regions of the world, recent studies showed that NDVI is effective in monitoring postfire vegetation succession [43,44]; however, our results indicate that NBR is the most reliable index for interannual comparisons and determining changes in the phenological pattern in this region of Bolivia, which allows for the detection of postfire regeneration. Fire severity levels based on dNBR and CBI indices are reliable methodologies that allow for determining the severity and dynamics of changes in postfire regeneration levels in forested and nonforested areas, as well as in their transitional stages. In the case of the dNBR, it is necessary to consider some important aspects, such as comparisons between the same stations. In Chiquitania, fires generally occur during the dry season [45,46]; after the first rains, there is a higher concentration of moisture in the soil, and vegetation resprouts develop, which could be interpreted as an apparent complete regeneration. In this study, we compared differences in the severity levels of dNBR after 1-2 years, considering seasonality. The dNBR index decreases over time as an effect of the disappearance of fire disturbance. In the case of CBI, the size of the plots and the practicality in data collection render it a replicable methodology to any site in the Chiquitania. Furthermore, we found high and significant correlation between the CBI and dNBR indices for 2019. The majority of studies showed medium-to-high correlation between CBI/GeoCBI severity levels and NBR-dNBR values [47][48][49][50][51]. However, all these correlation analyses between the two indices were performed in boreal ecosystems, and information for tropical ecosystems is almost nonexistent, so our research contributes to determining the reliability of using severity analysis based on remote sensing. In addition, our research was based on the classification proposed by Key and Benson [28], who defined the thresholds on the basis of research in boreal ecosystems, so that future research should be directed to define the categories of forest fire severity in tropical ecosystems.
Differences in fire severity levels in the same area are not fully understood. Some aspects that should be evaluated in future studies to assess these differences should consider the type of fire (ground, surface or crown fires), meteorological conditions, the moisture concentration of the biomass, and the type of forest fuel [52]. In certain forested regions of the Chiquitania, some fire-adapted species (e.g., bamboo Guadua paniculata), generates high levels of fine biomass, and increases fire intensity and forest flammability [53,54]. In the nonforested areas adapted to fire in the Laguna Marfil and Ñembi Guasu (e.g., Cerrado, Abayoy, natural savannah), the fires consumed the greatest amount of fine and medium biomass; in the forested areas of the four study areas, they were surface fires, but in certain areas of the forest, the appearance of crown fires was evident, demonstrating higher levels of severity dNBR and CBI, but without reaching the most extreme values recorded. In Ñembi Guasu, the protected area where a series of magnitude fire events with a high magnitude were recorded [55,56], we determined that severity levels were the most extreme for 2019, due to the high intensity of wildfires [3].
A fundamental topic in spatial forest ecology that is related to fire is mapping and modeling the complexity of postburn forest patterns and their changes over time [20]. The results yielded in our research show that the forested and nonforested ecosystems of the Chiquitania impacted by fires, are in the process of natural regeneration, whereas the only human interventions were the avoidance of degradation processes. Therefore, we recommend a more extensive protection and monitoring of the affected areas to ensure a successful natural regeneration and to avoid a recurrence of fires, which would reduce a successful regeneration of the areas impacted by the 2019-2020 fires. Although the results of the present work do not yet allow f or us to determine the time of complete recovery of the forests in terms of composition and structure because the growth of the tree canopy occurs much more slowly than the recovery in vegetation indices does. For this reason, it is necessary to continue monitoring the sites on the basis of remote-sensing and long-term field evaluations to better understand the dynamics, trends of change over time, and the ability of the different ecosystems of the Chiquitania to recover.

Conclusions
After the 2019 and 2020 fires, a process of natural regeneration of the ecosystems is evident among the four study sites. However, some regions are in a state of degradation, mainly in areas where the severity of the fires was high (Ñembi Guasu) and where these fire events were repeated in two consecutive years (Laguna Marfil and Copaibo). These results should be complemented with further studies, since it is necessary to continue developing long-term monitoring mechanisms that allow for us to deepen our understanding of fire ecology and natural regeneration processes in these sites. Our results demonstrate the importance of using monitoring technologies including remote sensing assessments in combination with targeted and cost-effective field work to obtain more accurate information on the responses of impacted ecosystems, both in terms of their ecological characteristics and the degree of intensity and repetitiveness of wildfires.
We recommend precautionary restoration action for areas affected by fires, and the establishment of safeguards and guarantees for the development of natural regeneration processes, as an action properly delimited, planned, and implemented according to the respective legal, technical, and administrative guidelines. Lastly, there is a clear need to establish a scientific and technical research agenda related to methodologies and analyses to evaluate scenarios in the framework of restoration.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/fire5030070/s1, Figure S1: Monthly phenological patterns of burned and unburned areas based on Sentinel-2 NDVI for 2019 and 2020 for the four study areas. Table S1: Fire severity levels in hectares and percentages (in parentheses) based on Sentinel-2 dNBR for 2019 and 2020 for the four study areas; Table S2: Fire severity levels based on CBI plots in the four study areas of the department of Santa Cruz.