*Article* **Remote Sensing Calculation of the Influence of Wildfire on Erosion in High Mountain Areas**

**Yolanda Sánchez Sánchez 1,\* , Antonio Martínez Graña <sup>1</sup> and Fernando Santos- Francés 2**


**Abstract:** Soil erosion is one of the most important environmental problems of the moment, especially in areas affected by wildfires. In this paper, we study pre-fire and post-fire erosion using remote sensing techniques with Sentinel-2 satellite images and LiDAR. The Normalized Burn Ratio is used to determine the areas affected by the fire that occurred on 18 August 2016 in the Natural Reserve of Garganta de los Infiernos (Cáceres). To calculate the erosion, the multi-criteria analysis is carried out from the RUSLE. Once all calculations were performed, there was a considerable increase in sediment production from 16 June 2016 (pre-fire) with an erosion of 31 T/ha·year to 16 June 2017 of 74 T/ha·year for areas of moderate fire severity, and an increase from 11 T/ha·year in 2016 to 70 T/ha·year for areas with a very high severity. From the NDVI, it was possible to verify that this also affected the recovery of post-fire vegetation, decreasing the NDVI index 0.36 in areas of moderate severity and 0.53 in areas of very high severity.

**Keywords:** vegetation dynamics; RUSLE; sentinel-2; soil erosion; wildfire

## **1. Introduction**

Research on soil erosion has been carried out for decades [1–3] because it is perceived as one of the most important environmental problems in the world, especially in high mountain regions with a Mediterranean climate since they have a high rainfall intensity, frequent outcrops of soft and weatherable rocks and scarce vegetation cover. This combination constitutes a favorable framework for natural soil erosion, which is aggravated by strong human pressure from certain activities that imply a change in vegetation, land use and soil properties. Consequently, soil loss increases exponentially giving rise to so-called anthropogenic erosion [4]. Given this problem, a particularly pressing issue is the quantification of erosion; different methods have been studied for its calculation: MUSLE [5], PESERA [6,7], Eurosem [8], Artificial Neural Network [9] and USLE [10]. In recent years, these methods have been integrated into Geographic Information Systems (GIS) to increase the accuracy of the calculation [11].

One of the most important causes in the acceleration of erosion rates are forest fires [12] since they represent a very abrupt disappearance of the vegetation cover, leaving the soil bare for weeks or months, increasing the production of runoff and sediment [13]. For this purpose, different indexes have been studied to evaluate the damage caused by forest fires, among them the NBR (Normalized Burn Ratio) [14], applied to post-fire mapping [15], from low resolution satellites [16].

For the calculation of erosion the most used method is the USLE (Universal Soil Loss Equation) or the RUSLE (Revised Universal Soil Loss Equation) with the disadvantage of using the latter for forest fires being that the dynamic calculation of the C factor is very complicated and very costly because it is necessary to revisit the study area continuously, and taking into account the canopy and understory vegetation, together with the Covered

**Citation:** Sánchez Sánchez, Y.; Martínez Graña, A.; Santos- Francés, F. Remote Sensing Calculation of the Influence of Wildfire on Erosion in High Mountain Areas. *Agronomy* **2021**, *11*, 1459. https://doi.org/ 10.3390/agronomy11081459

Academic Editor: Massimo Fagnano

Received: 23 June 2021 Accepted: 20 July 2021 Published: 22 July 2021

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

**Copyright:** © 2021 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/).

Cavity Fraction being very difficult to study with an accuracy of 10 m pixels. Therefore, different calculation methods have been used for this factor: vegetation cover [17], a radiometer [18] or from NDVI applied to other fields such as erosion in war zones [19]. different calculation methods have been used for this factor: vegetation cover [17], a radi‐ ometer [18] or from NDVI applied to other fields such as erosion in war zones [19]. Remote sensing has been used in other aspects related to wildfire: calculation of fuels

and taking into account the canopy and understory vegetation, together with the Covered Cavity Fraction being very difficult to study with an accuracy of 10 m pixels. Therefore,

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 2 of 13

Remote sensing has been used in other aspects related to wildfire: calculation of fuels models [20], post-fire recovery [21] and multitemporal study [22,23]. Satellite imagery has not only been used in wildfire but also in themes such as altimetry [24] bathymetry [25,26] and ice [27]. models [20], post‐fire recovery [21] and multitemporal study [22,23]. Satellite imagery has not only been used in wildfire but also in themes such as altimetry [24] bathymetry [25,26] and ice [27]. In this work, we related pre‐ and post‐forest fire erosion production from remote

In this work, we related pre- and post-forest fire erosion production from remote sensing, observing the areas with a remarkable increase in erosion and those more susceptible to sediment production. For this purpose, two images obtained from the Sentinel 2 Satellite with a time difference of 1 year (June 2016–June 2017) were compared from the data of vegetation vigorousness and land uses; combining it with LiDAR data and the Digital Terrain Model (DTM) allowed us to calculate with a high accuracy the erosion from the multivariate model of the RUSLE. The study also intended to differentiate the different severity of the Garganta de los Infiernos forest fire in September 2016, from calculations obtained by remote sensing and to see how this evolved over a period of time with respect to vegetation, erosion and sediment production. sensing, observing the areas with a remarkable increase in erosion and those more sus‐ ceptible to sediment production. For this purpose, two images obtained from the Sentinel 2 Satellite with a time difference of 1 year (June 2016–June 2017) were compared from the data of vegetation vigorousness and land uses; combining it with LiDAR data and the Digital Terrain Model (DTM) allowed us to calculate with a high accuracy the erosion from the multivariate model of the RUSLE. The study also intended to differentiate the different severity of the Garganta de los Infiernos forest fire in September 2016, from cal‐ culations obtained by remote sensing and to see how this evolved over a period of time with respect to vegetation, erosion and sediment production.

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

#### *2.1. Study Area 2.1. Study Area*

The study area is the Garganta de los Infiernos nature reserve in the Jerte Valley in the north of the province of Cáceres (Spain), with an area of 73 km<sup>2</sup> . It is predominantly steep slopes, starting from an altitude of 526 m to 2337 m. The main permanent watercourse is the stream of the "garganta de los infiernos". The average annual rainfall is 1300 mm, mainly between October and June. The lithology consists mainly of poorly developed soils due to the high slopes and the large influx of granite rocky outcrops. The study area is the Garganta de los Infiernos nature reserve in the Jerte Valley in the north of the province of Cáceres (Spain), with an area of 73 km2. It is predominantly steep slopes, starting from an altitude of 526 m to 2337 m. The main permanent water‐ course is the stream of the "garganta de los infiernos". The average annual rainfall is 1300 mm, mainly between October and June. The lithology consists mainly of poorly developed soils due to the high slopes and the large influx of granite rocky outcrops.

The forest fire (Figure 1) studied took place on 18 August 2016 and remained at level 2 until 20 August; it was extinguished on 28 August having finally razed approximately 1000 ha. The fire occurred in a mountainous area with a large influx of granitic rock formations; the burned vegetation was scrubland and rocky vegetation. The forest fire (Figure 1) studied took place on 18 August 2016 and remained at level 2 until 20 August; it was extinguished on 28 August having finally razed approximately 1000 ha. The fire occurred in a mountainous area with a large influx of granitic rock for‐ mations; the burned vegetation was scrubland and rocky vegetation.

**Figure 1.** Situation map of Garganta de los Infiernos (Study area) of Digital Terrain Model (DTM) with the perimeter of forest fire.

#### *2.2. Method 2.2. Method*

with the perimeter of forest fire.

The methodology followed is summarized in Figure 2 and will be explained step by step in the following subsections. The methodology followed is summarized in Figure 2 and will be explained step by step in the following subsections.

**Figure 1.** Situation map of Garganta de los Infiernos (Study area) of Digital Terrain Model (DTM)

*Agronomy* **2021**, *11*, x FOR PEER REVIEW 3 of 13

#### **Figure 2.** Workflow.

#### **Figure 2.** Workflow. *2.3. Remote Sensing*

2.3.1. Remotely Sensed Image from Sentinel 2

*2.3. Remote Sensing* 2.3.1. Remotely Sensed Image from Sentinel 2 Sentinel‐2 is a mission of the European Space Agency composed of the Launch of two satellites: Sentinel 2A (launched in June 2015) and Sentinel 2B (Launched in March 2017). This satellite was chosen because it is the European satellite providing the highest spatial resolution (10 m) among those offering free services [28] and a temporal resolution of five Sentinel-2 is a mission of the European Space Agency composed of the Launch of two satellites: Sentinel 2A (launched in June 2015) and Sentinel 2B (Launched in March 2017). This satellite was chosen because it is the European satellite providing the highest spatial resolution (10 m) among those offering free services [28] and a temporal resolution of five days. The satellite has 13 spectral bands (Table 1) ranging from visible and near-infrared (VNIR) to shortwave infrared (SWIR) wavelengths over an orbital swath of 290 km [29].

days. The satellite has 13 spectral bands (Table 1) ranging from visible and near‐infrared **Table 1.** Radiometric and Spatial resolution of Sentinel-2.


9 Water vapor 945 60 10 Cirrus detection 1375 60 11 SWIR 1610 20 12 SWIR 2190 20 In this work, we have studied 3 Sentinel-2 satellite images of the dates: 16 June 2016 (Pre-fire), 4 September 2016 (Post-fire) and 2 June 2017 (spring post-fire). These dates were chosen because they were the images with the best conditions for obtaining the necessary indices. In addition, during the months of May and June, field campaigns were carried out for the mapping of fuel models [20], so the data were used to corroborate the calculated indices. The post-fire image chosen was in September since it was the first available image

8a NIR 865 20

without smoke, without clouds and with high quality. So it is the closest available image to the date of interest.

From these 3 images, different band combinations were made between them: RGB Geology (band 12, 11 and 2) and RGB Agriculture (band 11, 8 and 2) which have been used for soil classification. Two indices were also calculated: NDVI (Normalized Difference Vegetation Index) [30] and NBR (Normalize Burn Ratio) [31].

#### 2.3.2. LiDAR

LiDAR (Light Detection and Ranging) [32] information from Spanish National Geographic Institute of the PNOA (National Aerial Orthophotography Plan) with a resolution of 0.5 points/m was used, and from this information, and by interpolation, the Digital Elevation Model with a resolution of 5 m was obtained. The data of the study area were obtained in 2010 and the images had a calibration of the LiDAR sensor, a maximum of five returns per pulse, and a pre-classification of said returns. The flight height of the sensor was a maximum of 3000 m from the ground with a horizontal accuracy of 0.30 m and vertical accuracy of 0.20 m.

#### *2.4. Normalized Burn Ratio (NBR)*

This index allows us to delimit the perimeter of the wildfire [33] and also to delimit the areas of the fire according to its severity with great precision [34]. Fire severity is a critical factor because of its direct relationship with the biomass consumed, so it is linked to post-fire vegetation and hydrogeomorphological recovery [14]. This index integrates the two bands that best show the combustion B8 (NIR-842 nm) and B12 (SWIR-2190 nm) in the following Equation (1)

$$\text{NBR} = (\text{NIR} - \text{SWIR}) / (\text{NIR} + \text{SWIR}) \tag{1}$$

This gives the NBR for a given date. In order to appreciate the change derived from forest fire (dNBR), the post-fire NBR is subtracted from the pre-fire NBR (Equation (2))

$$\text{NBRR} = \text{NBRPrefire} - \text{NBRpostfire} \tag{2}$$

Once this index is obtained, the study area is classified into different zones depending on the forest fire severity [35]. For this purpose, we will reclassify the output raster into 4 classes: high severity (dNBR > 680), moderate severity (dNBR > 275), low severity (dNBR > 90) and unburned (dNBR < 90) [16].

#### *2.5. RUSLE (Revised Universal Soil Loss Equation)*

In this work, the RUSLE (Revised Universal Soil Loss Equation) [36] was used to quantify the erosion rate of the Garganta de los Infiernos before the fire, just after the fire and one year after the fire. This equation was the most appropriate for this calculation since it is the one that has the widest application as it is a parametric and totally empirical model. It is a totally empirical formula (Equation (3)) that attempts to interpret the erosive mechanisms by their causes and effects.

$$\mathbf{A} = \mathbf{R} \times \mathbf{K} \times \mathbf{L} \mathbf{S} \times \mathbf{C} \times \mathbf{P} \tag{3}$$

where A is the Average annual soil loss (T/ha·year), R is the rainfall erosivity (MJ·mm/ha·h·year), K is the soil erodibility (Tm/ha h/MJ·mm), LS is the hill slope length and steepness (dimensionless), C is the vegetation factor (dimensionless) and P is the support practice (dimensionless).

• Rainfall erosivity factor (R)

Rainfall erosivity factor is the factor that represents the effect of rainfall on soil erosion. It is calculated as the result of multiplying the kinetic energy of rainfall by the maximum intensity during 30 min of precipitation [37] (Equation (4))

$$\mathbf{R} = \mathbf{E} \cdot \mathbf{I}\_{\text{30}} \tag{4}$$

where R is the rainfall erosivity (MJ·mm/ha·h·year), E is total storm energy (Mj/ha·year) and I<sup>30</sup> the maximum 30 min rainfall intensity (mm/hour).

These data were obtained from the rainfall intensity recorded every 30 min from the Aldehuela del Jerte (CC04) and Valdeastillas (CC17) stations and an extrapolation was performed to calculate the R value at the peak with the highest altitude in the study area (Table 2) for the month of June 2016, the month of September 2016 and the month of June 2017.

**Table 2.** Location of weather stations and data of the R factor (Information elaborated using that obtained from the AEMET (Agencia Estatal de Meteorología).


• Soil erodibility factor (K)

This factor represents the response of the soil to a given erosive force or mechanism, or, in other words, the susceptibility of the soil to erosion [38]. It represents the quantified soil loss per unit of erosivity in a standard plot 22.6 m long with a slope of 9%. The data were obtained from the experimental plots of the national soil inventory of Cáceres [39] and an unsupervised classification of the RGBGeology and RGBAgriculture band composition was performed to separate the zone into rocky, bare soil and vegetation. The rocky zone was assigned the value 0.134 (mean measured R value for the plutonic rock lithofacies). For the other zones, the R value of the experimental plots was extrapolated.

• Topographic factor (LS)

The LS factor responds to the combined effect of slope length and slope angle. Its value is used to estimate the soil losses that occur on a sloping terrain compared to the losses per unit if the same rainfall were to fall on a standard plot with identical soil type, crop and management conditions.

This factor was calculated from the DTM-LiDAR with a resolution of 5 m; from the DTM, the slope and the flow direction were calculated to subsequently calculate the flow accumulation, and with the following formula the LS factor was calculated (Equation (5)) [40]

LS = (flow accumulation cell size/22.13)0.4 (sin ([slope] <sup>×</sup> 0.01745)/(0.0896)1.3 (5)

• Vegetation factor (C)

Factor C takes into account vegetation cover to protect the soil. Vegetation is a relevant factor in erosion, vegetation cover reduces the energy of rainfall by intercepting it and prevents it from falling directly to the ground in addition to promoting the infiltration of runoff water [41]. For this article, the C factor was calculated from the NDVI. This index is effective for quantifying green vegetation; it normalizes the scattering of green leaves in the near-infrared wavelength and the absorption of chlorophyll in the red wavelength [30]. The index has been used to see the fluctuations of vegetation at different times of the year and the difference of pre- and post-fire vegetation [18] by applying the formula (Equation (6)) [2]

$$\mathcal{C}\_{\text{factor}} = \exp[-\alpha \,\,\frac{\text{NDVI}}{\text{ $\mathfrak{F}$ } - \text{NDVI}}] \tag{6}$$

where α and β parameters determine the shape of the NDVI curve. Reasonable results are produced using values of α = 2 and β = 1.

• Conservation practices (P)

Conservation practices are measures to reduce runoff, generally carried out in cultivated areas. In the study area, there are no cultivated areas, so conservation practices are discouraged.

#### *2.6. Validation*

Finally, 200 sampling points randomly distributed over the entire surface of the study area were used, from which the C-factor data obtained from the NDVI index by remote sensing were obtained and compared with the actual C-factor data and with the data obtained in the field. The mean square error (MSE) index was used to verify the validation of the cartography (Equation (7)).

$$\text{MSE} = \frac{1}{n} \sum\_{i=1}^{n} (\text{\textdegree c} - \text{\textdegree r})^2 \tag{7}$$

where *n* is the number of data points, Xc the value of the generated cartography and Xr the actual value for the data point.

#### *2.7. Statistical Analysis*

The correlation methods were used because they are the statistical tools that indicate the relationship between two variables. In this case, the variables that had a more direct affect in post-fire erosion were studied. The correlation method used was the Pearson coefficient that measures the degree of covariance between different linearly related variables. Therefore, a linear correlation procedure was performed to examine the relationship between the covariance of throat erosion post-fire with NDVI, factor C and factor R and to determine the relationship of each variable. The SPSS Statistics 25 software was used to carry out these analyzes.

#### **3. Results**

The process carried out to calculate the relationship between the pre-fire erosion of the Garganta de los Infiernos and the post-fire erosion was based on remote sensing data, so that different erosion values were obtained depending on the different fire severity zones.

In the first step of the process, calculation of the severity of the fire (Figure 3), four zones could be distinguished: where the fire had not affected (80%), zones near or with occurrence of other minor fires (7%) with a lower severity, zones of the perimeter where the severity of the fire has been moderate (10%) and the central focus of the fire with a high severity (3%).

Once the fire severity zones boundaries were marked, the NDVI index was calculated for the different dates (16 June 2016, 5 September 2016 and 2 June 2017) (Figure 4a–c) where the influence of the fire on the vegetation and subsequent regrowth of the vegetation was clearly verified.

**Figure 3.** Map of areas of forest fires severity by the NBR and occupancy graph by severity. **Figure 3.** Map of areas of forest fires severity by the NBR and occupancy graph by severity.

Once the fire severity zones boundaries were marked, the NDVI index was calculated for the different dates (16 June 2016, 5 September 2016 and 2 June 2017) (Figure 4a–c) **Table 3.** Table of the output of the SPSS with the results of the Pearson correlation between the factors that cause erosion and soil loss.


the fire severity had been very high: 0.55 pre‐fire and 0.03 post‐fire. Vegetation recovery also depended on the severity of the fire. In the medium severity zone, the difference in \*\* The correlation is significant at the level 0.01 (bilateral); The correlation is significant at the level 0.05 (bilateral).

NDVI from 2016 to 2017 was 0.36, while in the high severity zone it was 0.52, which shows that vegetation recovery was slower in this zone. Once the NDVI values were obtained, the C‐factor mapping was created, this map‐ ping was compared with the C‐factor mapping produced from the values established by the Soil Conservation Service of the United States and with the values of the Second Na‐ tional Forest Inventory for the province of Cáceres [39]. The comparison of both C‐factor mappings gave an MSE of 0.15, which indicates that the procedure for obtaining the C‐ In Figure 4d, it can be seen that in the area where the fire had no affect, the NDVI did not vary, however, in the area where the influence of the fire was low the NDVI varied from 0.52 pre-fire to 0.33 post-fire. Where the change is most noticeable is in the area where the fire severity had been very high: 0.55 pre-fire and 0.03 post-fire. Vegetation recovery also depended on the severity of the fire. In the medium severity zone, the difference in NDVI from 2016 to 2017 was 0.36, while in the high severity zone it was 0.52, which shows that vegetation recovery was slower in this zone.

factor mapping from the satellite images was highly accurate. From the mapping of the C factor, R factor, LS factor and K factor, the RUSLE was calculated and which shows a clear difference between the month of June 2016 and the month of June 2017, in which the mapping of orography and lithology were the same; the mapping of climatology varied since the month of June 2017 was a rainier month than the month of June 2016. However, the most pronounced change was the vegetation, which, in Once the NDVI values were obtained, the C-factor mapping was created, this mapping was compared with the C-factor mapping produced from the values established by the Soil Conservation Service of the United States and with the values of the Second National Forest Inventory for the province of Cáceres [39]. The comparison of both C-factor mappings gave an MSE of 0.15, which indicates that the procedure for obtaining the C-factor mapping from the satellite images was highly accurate.

the areas of high severity and moderate severity of the fire, almost disappeared. Therefore, a comparison was made between the erosion mapping before the fire and after the fire. In Figure 5 (Figure 5a,b), it can be seen that in the study area not affected by the fire, erosion was similar, while in the study area that was affected by the fire (Figure 5b,d) there are areas that went from low to high or very high erosion. The moderate severity zone in 2016 possessed an average erosion of 31 T/ha∙year, while in 2017 it possessed an average erosion of 74 T/ha∙year, 2.5 times higher. Even more remarkable was the change From the mapping of the C factor, R factor, LS factor and K factor, the RUSLE was calculated and which shows a clear difference between the month of June 2016 and the month of June 2017, in which the mapping of orography and lithology were the same; the mapping of climatology varied since the month of June 2017 was a rainier month than the month of June 2016. However, the most pronounced change was the vegetation, which, in the areas of high severity and moderate severity of the fire, almost disappeared. Therefore, a comparison was made between the erosion mapping before the fire and after the fire.

**Figure 4.** (**A**) NDVI 16 June 2016; (**B**) NDVI 5 September 2016; (**C**) NDVI 2 June 2017; (**D**) Graph of Table 3. the covariation of vegetation and climatology factors versus erosion data in 2017 was studied. The results of Pearson correlation between C factor, NDVI, R factor, erosion and fire severity are presented in Table 2. **Figure 4.** (**A**) NDVI 16 June 2016; (**B**) NDVI 5 September 2016; (**C**) NDVI 2 June 2017; (**D**) Graph of Table 3. the covariation of vegetation and climatology factors versus erosion data in 2017 was studied. The results of Pearson correlation between C factor, NDVI, R factor, erosion and fire severity are presented in Table 2.

The analysis shows that there was an inverse correlation between the C factor and erosion with a high significance, the same correlation value gives us the relationship be‐ tween NDVI and erosion since the C factor was obtained from NDVI, which indicates a precision in the calculations. The R factor and RUSLE had a positive covariance, which indicates that the relationship between rainfall intensity and erosion was positive, but the significance, despite being high, was lower than in the case of vegetation; however, in areas where the soil was not affected by forest fire, rainfall was more important than veg‐ The analysis shows that there was an inverse correlation between the C factor and erosion with a high significance, the same correlation value gives us the relationship between NDVI and erosion since the C factor was obtained from NDVI, which indicates a precision in the calculations. The R factor and RUSLE had a positive covariance, which indicates that the relationship between rainfall intensity and erosion was positive, but the significance, despite being high, was lower than in the case of vegetation; however, in areas where the soil was not affected by forest fire, rainfall was more important than vegetation.

in the areas of high fire severity: in 2016 these areas had an erosion of 11 T/ha∙year, how‐

ever, in 2017, they reached an erosion of 70 T/ha∙year, more than six times higher.

etation. **4. Discussion** The NBR index has been used in several forest fire studies with lower resolution sat‐ ellites [42], or with Sentinel‐2 [43]. In all of them, the relationship between fire severity [44] and remotely sensed data from this index has been very good. The resulting difference between the 2016 and 2017 NDVI is due to the fact that, after In Figure 5 (Figure 5a,b), it can be seen that in the study area not affected by the fire, erosion was similar, while in the study area that was affected by the fire (Figure 5b,d) there are areas that went from low to high or very high erosion. The moderate severity zone in 2016 possessed an average erosion of 31 T/ha·year, while in 2017 it possessed an average erosion of 74 T/ha·year, 2.5 times higher. Even more remarkable was the change in the areas of high fire severity: in 2016 these areas had an erosion of 11 T/ha·year, however, in 2017, they reached an erosion of 70 T/ha·year, more than six times higher.

a wildfire, vegetation biomass and wattle layer are converted into ash, carbon and organic matter altered by the fire, which are partially carried away by runoff [45] and partially incorporated into the soil [46], retarding the growth of vegetation cover by modifying re‐ growth dynamics. The vegetation did not fully recover one year afterthe fire, thus limiting

the input of waffles [47].

and that soil water repellency.

This analysis shows us that the increased erosion rate is due to the changes that wildfire ignition causes in soils and vegetation. Different authors have shown that the increase in temperatures up to 550 °C from a wildfire completely destroys soil hydrophobicity [48]

**Figure 5.** (**A**) Erosion map 16 June 2016; (**B**) Zoom of erosion map in forest fire area 16 June 2016; (**C**) Erosion map 02 June 2017; (**D**) Zoom of erosion map in forest fire area 02 June 2017; (**E**) Comparative graph of the erosion in the different months by area of forest fire severity. **Figure 5.** (**A**) Erosion map 16 June 2016; (**B**) Zoom of erosion map in forest fire area 16 June 2016; (**C**) Erosion map 02 June 2017; (**D**) Zoom of erosion map in forest fire area 02 June 2017; (**E**) Comparative graph of the erosion in the different months by area of forest fire severity.

#### **4. Discussion**

The NBR index has been used in several forest fire studies with lower resolution satellites [42], or with Sentinel-2 [43]. In all of them, the relationship between fire severity [44] and remotely sensed data from this index has been very good.

The resulting difference between the 2016 and 2017 NDVI is due to the fact that, after a wildfire, vegetation biomass and wattle layer are converted into ash, carbon and organic matter altered by the fire, which are partially carried away by runoff [45] and partially incorporated into the soil [46], retarding the growth of vegetation cover by modifying regrowth dynamics. The vegetation did not fully recover one year after the fire, thus limiting the input of waffles [47].

This analysis shows us that the increased erosion rate is due to the changes that wildfire ignition causes in soils and vegetation. Different authors have shown that the increase in temperatures up to 550 ◦C from a wildfire completely destroys soil hydrophobicity [48] and that soil water repellency.

Decreases enhancing splash erosion processes [49]. If to that is added the direct relationship with rainfall intensity, then this leads to an increase in runoff rates reducing the availability of nutrients such as organic carbon and water [50]. In other similar studies with similar vegetation and climate conditions, they concluded [7,43], as in this article, that the erosion rate increases greatly during the first two years after the fire and that the vegetation would largely recover after the fifth year [51]. Thus, we can highlight an exponential increase in erosion in mountain areas after the occurrence of a fire due to the loss of vegetation and the exposure of bare soil in possible high intensity rainfall events, and that, even after the following period of maximum vegetation reproduction, erosion is very high [52].

#### **5. Conclusions**

Once the results were analyzed, it wase concluded that soil loss for the same area, with the same lithology and orography, is much more abundant after a fire than after a burn, which implies soil and vegetation degradation. There is a clear relationship between the occurrence of a fire and the increase in soil loss, which can be up to 6.5 times higher in areas with a high fire affectation, as a consequence of the massive loss of vegetation. This could lead to contamination of downstream aquifers due to ash production, to saturation of reservoirs due to sedimentation and to a loss of soil productivity, which is why the subsequent treatment of a forest fire by the Spanish Ministry of the Environment is of great concern.

On the contrary, the NBR index is found to be effective in quantifying fire damage in the different affected areas and shows a very high spatial coincidence between the areas of highest fire severity and the areas of highest post-fire sediment production. It is also concluded that the NDVI index is very effective in calculating the RUSLE vegetation factor because it is able to calculate vegetation vigor dynamically, instantaneously and parametrically with high accuracy (MSE = 0.15). Using the NDVI index together with the NBR it can be seen that, in the areas most impacted by the fire, the vegetation takes longer to recover since the damage was greater, so the erosion rate increases linearly at the same time that the vegetation has decreased and the intensity of rainfall increases.

Finally, with remote sensing from Sentinel-2 images, it is possible to quantify and compare the erosion of a pre- and post-fire area in the different areas of severity of forest fires with high accuracy and at low cost. This methodology allows us to make a useful tool to map a burned area with a high spatial resolution (10 × 10 m) and to study the evolution of post-fire vegetation with a temporal resolution of only 10 days, which allows us to make better decisions when carrying out the soil study.

**Author Contributions:** Conceptualization, Y.S.S. and A.M.G.; methodology, Y.S.S.; software, Y.S.S.; validation, Y.S.S. and F.S.-F.; formal analysis, A.M.G. and F.S.-F.; investigation, Y.S.S. and F.S.-F.; resources, Y.S.S. and A.M.G.; writing—original draft preparation, Y.S.S.; writing—review and editing, A.M.G.; supervision, F.S.-F.; project administration, A.M.G.; funding acquisition, A.M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the project SA044G18 of Regional Government of Castilla y Leon, and the GEAPAGE research group (Environmental Geomorphology) of the University of Salamanca.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research was funded by the project SA044G18 of Regional Government of Castilla y Leon, and the GEAPAGE research group (Environmental Geomorphology) of the University of Salamanca.

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

## **References**

