Combining Methods to Estimate Post-Fire Soil Erosion Using Remote Sensing Data

: The increasing number of wildﬁres in southern Europe is making our ecosystem more vulnerable to water erosion; i.e., the loss of vegetation and subsequent runoff increase cause a shift in large quantities of sediment. Fire severity has been recognized as one of the most important parameters controlling the magnitude of post-ﬁre soil erosion. In this paper, we adopted a combination of methods to easily assess post-ﬁre erosion and prevent potential risk in subsequent rain events. The model presented is structured into three modules that were implemented in a GIS environment. The ﬁrst module estimates ﬁre severity with the Monitoring Trends in Burn Severity (MTBS) method; the second estimates runoff with rainfall depth–duration curves and the Soil Conservation Service Curve Number (SCS-CN) method; and the third estimates pre- and post-ﬁre soil erosion. In addition, two post-ﬁre scenarios were analyzed to assess the inﬂuence of ﬁre severity on soil erosion: the former based on the Normalized Difference Vegetation Index (NDVI) and the latter on the Relative differenced Normalized Burn Index (RdNBR). The results obtained in both scenarios are quite similar and demonstrate that transitional areas, such as rangelands and rangelands with bush, are the most vulnerable because they show a signiﬁcant increase in erosion following a ﬁre event. The study ﬁndings are of secondary importance to the combined approach devised because the focal point of the study is to create the basis for a future tool to facilitate decision making in landscape management.


Introduction
In the last decade, an increasing number of fires have occurred in southern Europe, causing severe damage to property, infrastructure, and natural resources [1]. The combination of climate change and human activities incrementing fire events has made our ecosystems more vulnerable to wildfires [2]. Post-fire soil erosion represents one of the most important detrimental indirect effects of fires in the Mediterranean countries of southern Europe. It has been well demonstrated that soil erosion and loss of slope stability increase as a result of fire occurrence [3][4][5][6][7].
Wildfires alter vegetation cover such as shrubs and trees, often causing severe damage to stems and total loss of canopy [2,8,9]. This may generate long-term impacts on the landscape because plants and trees lose their ability to stabilize the soil [10,11]. Root systems anchor the soil making it stable beneath the surface layer [12]. Low plants provide ground cover from wind, while higher trees limit the strength of rain before hitting the ground. Once plants are established, their life cycles help return nutrients to the soil to encourage future plant growth and maintain adequate moisture levels thus reducing soil erosion [13,14]. More specifically, after fire occurrence, there is a loss of nutrient-storage capacity in the

Study Area
The Apulia region represents the easternmost region in Italy and is located at a latitude of 39°50′-41°50′N and a longitude of 15°50′-18°50′E (Figure 1a). The climate is typically Mediterranean with mild rainy winters, dry hot summers and rainfall in the form of short but intense storms. It is precisely in summer that most of the wildfires occur, favored by these climatic conditions and moderate-to-high mean wind speed.   Apulia's landscape is characterized by a great variety of species and boasts a significant geological and natural heritage. Despite this variety, according to the National Forest Inventory (2015) [46], the forests in Apulia cover an area of 189,086 ha, which represents 9.7% of the territory, an amount significantly lower than the national average. About one sixth of this area was affected by fires in the period between 2005 and 2016 [47].
Our combining approach was tested on four study sites (Figure 1b) derived from the 2019 summer wildfire database (provided by the Apulia Region Civil Protection Department). The sites were selected on the basis of their slope characteristics and the presence of natural areas such as deciduous wood, transitional woodlands, natural grassland, bushes and shrubs, and sparsely vegetated areas. One site is located on the Tavoliere Plain, a large expanse covering 3000 km 2 and about 15% of the total area of the Apulia region ( Figure 1a). Geologically, it is part of a sedimentary basin originating in the early Pliocene during the final phases of orogenesis of the southern Apennines. Therefore, the outcropping soils essentially consist of marine sediments (Plio-Pleistocene clays) and alluvial deposits (gravels, sands, and silts) [48,49]. The study site can be considered a wildland urban interface, precisely located in the municipality of Lucera (LUC, Figure 1c). The LUC site is characterized by deciduous woods, transitional woodlands, natural grassland, and sparsely vegetated areas. The other three study sites are located on the Gargano promontory. This landscape is characterized by a relevant and heterogeneous geodiversity. It consists of extrusions of carbonate rocks in different sedimentary contexts, strongly modeled by tectonic systems [50]. Gargano (Figure 1a) is also considered a hot spot in terms of wildfire occurrence and frequency [47,51,52]. The study sites, located in the municipalities of Carpino (CAR, Figure 1d) and Ischitella (ISC_1, Figure 1e and ISC_2, Figure 1f), are characterized by deciduous woods, transitional woodlands, natural grassland, pastures, bushes and shrubs, and olive groves.
The four study sites were burned by wildfires between July and August of 2019. The largest (39.5 ha) occurred in Lucera and affected a hillslope with slopes between 5% and 50%. The other three wildfire events affected smaller areas with slope not exceeding 30% (Table 1).

Model Description
The combining approach consists of three modules, based on the structure of the POSTFIRE model proposed by Fox et al. [40], which are:
Erosion-pre-and post-fire scenario

Fire Severity Classification Module
In this module, the methodology of the Monitoring Trends in Burn Severity (MTBS) program was used for the identification of the burned areas and the classification of fire severity. This method provides a classification that considers the ecological significance and the site dependence of fire severity. As a consequence, the thresholds of fire severity classes are defined for each study site following the instructions provided by the MTBS program [53][54][55]. According to the aim of our study, an initial assessment was performed acquiring images within a short time pre-and post-fire event [54,56]. Firstly, the fire events were identified from the wildfire point database of the Apulia Region Civil Protection Department. Subsequently, the perimeters of the burned areas were identified starting from these punctual data and using different tools, as provided in the MTBS methodology. The analyses were performed using the land use map of Apulia Region derived from CORINE Land Cover [57], high-resolution True Color Image (TCI) from Google Earth and Esri World Imagery [58,59], and Level-1C products of a Sentinel-2 (S-2) satellite [60]. The S-2 images were cloud-free (Supplementary Table S1) and subjected to Top of the Atmosphere corrections. By combining different wavebands, the following indices and images were obtained: where: where NIR is the near-infrared band and SWIR is the shortwave infrared band [56]; • differenced Normalized Difference Vegetation Index (dNDVI): where: where NIR represents the near-infrared band and RED the red wavebands [61]; • TCI from Sentinel-2 Level-1C products; • pre-and post-fire SWIR composite images.
Lastly, fire severity was analyzed to identify the significant thresholds for each wildfire useful for discriminating the 4 severity classes: unburned, low, moderate, and high, as provided in the MTBS methodology. In addition to the products previously listed, the Relative differenced Normalized Burn Index (RdNBR) [62] was used for the interpretation of pixels subject to high severity: RdNBR = dNBR/(SquareRoot (ABS(NBR pre-fire ))) (5)

Runoff Module
In the second module, the Soil Conservation Service Curve Number (SCS-CN) method [63,64] was used for computing the depth of surface runoff (RO, mm) for a given rainfall event both under pre-and post-fire conditions. Consequently, the module was divided into two submodules: the first aimed at rainfall analysis and the second at calculating runoff.
The meteorological stations of interest for the study sites were identified by building Thiessen polygons [65] of the Apulia networks of the Regional Agency for Irrigation and Forestry Activities (ARIF) and Apulia Region Civil Protection Department. Time series of daily rainfall data were analyzed. After estimating the parameters using the moment method [66] and evaluating their fit with the Kolmogorov-Smirnov test [67], the data were modeled using the Gumbel distribution to obtain depth-duration-frequency curves [68][69][70]. The values of the depth of probable maximum precipitation (PMP) for a 1-day duration and 1.5-year return period [40] were extrapolated from the curves obtained.
According to the SCS-CN method, RO was evaluated using the following equations: where P (mm) is the value of daily rainfall with a 1.5-year return period; S (mm) is the potential infiltration; CN (dimensionless) is the curve number that ranges from 1 to 100 depending on the hydrological soil group and cover type, treatment, and hydrological condition.
Land use was harmonized and classified into the following macro-categories: deciduous forest, olive grove, rangeland, and rangeland with bush. The specific land use and area (percentage) for each case study are reported in Supplementary Table S2. In the pre-fire scenario, CN values provided by the Apulia section of the Southern Apennine District Basin Authority [71] were used ( Table 2). CNs were adjusted with the Antecedent Moisture Condition parameter (AMC) according to Hawkins et al. [72]. The AMC parameter depends on the total of 5 days preceding the rainfall event (three classes: I, II, and III representative of dry, average, and wet soil condition) and the season category (dormant or growing). The worst hydrological conditions occur when soil is almost saturated, and the excess precipitation is transformed into runoff [73]. For this reason, the AMCIII (highest soil moisture) class was considered (Supplementary Table S3). A-high, B-medium-high, C-medium-low, D-low permeability.
In the post-fire scenario, soil permeability decreased [74]. Based on much scientific literature, which uses simulation of different scenarios coupled with a land use update model to assess wildfire effects, we increased CNIIIs by 15, 10, 5, and 0 for the high, moderate, low, and unburned fire severity classes, respectively [42,75].

Erosion Module
The soil erosion module is based on the POSTFIRE model proposed by Fox et al. [40], which provides a combination of different previously created layers in a GIS environment to obtain both erosion pre-fire and post-fire maps. In the combining approach, to introduce the influence of fire severity on soil erosion, two different post-fire scenarios were proposed: NDVI-scenario and RdNBR-scenario.
Due to the absence of data measured directly in field, the sediment concentration values for the main land uses detected in the study areas were retrieved from the Soil and Water Assessment Tool (SWAT) [76]. The SWAT is a semi-distributed hydrological model that can be used to compute the streamflow and sediment transport in a watershed and to assess the results at different spatial and temporal scales [77]. Abdelwahab et al. [78] calibrated and validated the SWAT model in the Dauno Sub-Apennine area using continuously measured data of streamflow and sediment load. Moreover, the spatialized model results, in terms of sediment yield, were cross-validated by comparing them with the soil erosion data obtained by Panagos et al. [79] using the European model RUSLE2015 showing good agreement [78]. Hence, in this work the SWAT model sediment concentration spatialized findings were first subdivided into the different land uses and then into the main slope classes according to Fox et al. [40] (Table 3). Pre-fire and post-fire erosion values were obtained by combining runoff, slope classes and sediment concentration GIS layers. For each studied site, two post-fire scenarios were implemented. In the NDVI-scenario, an NDVI severity factor (NDVI SF ) was considered to correlate the change in vegetation cover with soil erosion [41,80]. This factor was obtained using the following Equation (8): The higher values are representative of the pixels where the change in vegetation cover in the post-fire scenario was higher.
In the second scenario, a RdNBR severity factor (RdNBR SF ) was considered to correlate fire severity with soil erosion. In this case initially, all the RdNBR negative values were set to 0 because they relate to an increase in vegetation cover and therefore do not contribute to an increase in erosion in the post-fire scenario [62]. Subsequently, a min-max normalization was applied to the RdNBR values and the RdNBR factor was obtained using the following Equation (9).

Fire Severity Classification
The fire severity classification ( Figure 2) showed a predominance of the moderate severity class (~30%-40% of the case study areas). The most fire-prone area is the LUC site where approximately 75% of the deciduous forest (>60% of the total LUC area) was affected by a low-to high-severity fire ( Figure 3).     Figure 3 shows that regardless of the area occupied in the site by rangelands (33% LUC, 46% CAR, 37% ISC_1, 11% ISC_2), this land use was the most affected by fire (80% LUC, 85% CAR, 95% ISC_1, 90% ISC_2). On the contrary, olive groves showed lower burnt area rates (35% ISC_1, 68% ISC_2), even when the land use is widespread as in the case study (63% ISC_1, 45% ISC_2).

Rainfall and Runoff
In the second module the PMP depth for a 1-day duration and 1.5-year return period was calculated for each study site (Table 4). Starting from these data and the application of the SCS-CN method, as previously described, the pre-and post-fire runoff values were obtained (Supplementary Figure S1). Variations in runoff values are due to both the fire severity class and the hydrological soil group. As shown in Table 5, the influence of these factors caused a variation between 10%-50% of runoff. The exceptions are olive groves, where the runoff variation was less than 10%, and two cases where runoff increased more than 50% (i.e., ISC_1-Rangeland, CAR-Deciduous forest).

Erosion
The results obtained in the third module are presented in Table 6. The average post-fire erosion values per land use in the pre-and post-fire scenarios (NDVI and RdNBR) within the study sites are reported. Higher values of average erosion, both pre-and post-fire, were observed in the rangelands. The differences between the average values of erosion in the NDVI and RdNBR scenarios were never more than 5% (order of hundredths), and the greatest values could be found in the rangelands. At the study site level, the greatest differences between the two scenarios were found in the site, as shown also in Figure 4. The figure shows the erosion values (t/ha) following a PMP for a 1-day duration and 1.5-year return period under a pre-fire condition and in the two post-fire scenarios (NDVI and RdNBR). In all study sites the results showed a tendency for eroded sediment to increase by 50% in the low-severity class, between 100% and 150% in the moderate class and between 150% and 300% in the high-severity class. As an example, Figure 5 illustrates the box plots of soil erosion values for each fire severity class and land use in the three scenarios of the LUC site.

Fire Severity Classification
Our study aim was to identify the areas most vulnerable to the action of atmospheric agents in the short term after fire events. For this reason, the first step of the model was the classification of fire severity, which is associated with the degree of immediate postfire environmental changes (biomass, fire products, and soil exposure) caused by fire and evaluated with an initial assessment [56,81,82]. In addition to supporting the initial assessment, the MTBS method allows the evaluation of fire effects on heterogeneous areas [54,55], such as those of the case studies ( Figure 3). The analysis of severity classification shows that the areal percentage of fire severity class per land use is site-dependent, i.e., amount of land use present in the site, border, or central location. However, rangelands demonstrate fire-prone behavior. As previously shown, rangelands are the land use most affected by fire. On the contrary, olive groves show lower burnt area rates even when this land use is widespread in the case study. These results are supported by D'Este et al. [83].
The MTBS model is robust, and the program has developed a tool which uses threshold ranges to identify fire severity classes for ecological systems that differ from those present in the Apulia region. Fire severity classification is, in fact, site-dependent. The identification of class thresholds is subjective [54], and the results obtained from the spectral indices are influenced by the vegetation characteristics as well [84]. Therefore, further studies on the severity of the green areas of Apulia are needed to define which spectral indices are most suitable for assessing severity and to obtain threshold ranges of severity classes. These are operations that would expedite severity classification.

Fire Severity Classification
Our study aim was to identify the areas most vulnerable to the action of atmospheric agents in the short term after fire events. For this reason, the first step of the model was the classification of fire severity, which is associated with the degree of immediate post-fire environmental changes (biomass, fire products, and soil exposure) caused by fire and evaluated with an initial assessment [56,81,82]. In addition to supporting the initial assessment, the MTBS method allows the evaluation of fire effects on heterogeneous areas [54,55], such as those of the case studies ( Figure 3). The analysis of severity classification shows that the areal percentage of fire severity class per land use is site-dependent, i.e., amount of land use present in the site, border, or central location. However, rangelands demonstrate fire-prone behavior. As previously shown, rangelands are the land use most affected by fire. On the contrary, olive groves show lower burnt area rates even when this land use is widespread in the case study. These results are supported by D'Este et al. [83].
The MTBS model is robust, and the program has developed a tool which uses threshold ranges to identify fire severity classes for ecological systems that differ from those present in the Apulia region. Fire severity classification is, in fact, site-dependent. The identification of class thresholds is subjective [54], and the results obtained from the spectral indices are influenced by the vegetation characteristics as well [84]. Therefore, further studies on the severity of the green areas of Apulia are needed to define which spectral indices are most suitable for assessing severity and to obtain threshold ranges of severity classes. These are operations that would expedite severity classification.

Rainfall and Runoff
The runoff module was structured to consider the worst conditions affecting post-fire erosion and therefore runoff. Determining rainfall thresholds for the runoff estimates is not sufficient given the many variables which influence the phenomenon. These variables partly depend on the variability of the site (e.g., burn severity, speed of vegetation recovery, sediment availability, and geomorphology) [85,86] and partly on rainfall. The literature has widely demonstrated that few (albeit intense) events favor major erosion phenomena [30,33,[87][88][89].
It also has been shown, however, that erosion depends on the space-time variability of the rains and on cumulative rainfall. In particular, intense rainfall combined with reduced infiltration leads to decreased hydraulic conductivity increasing runoff, especially in soils in which infiltration-excess mechanisms are prevalent, such as in the Mediterranean environment [90]. The increase in soil water repellence in the first period after fire occurrence is one of the possible causes of decreased infiltration capacity [91]. The AMC is also an aspect to consider. In fact, as stated by Todisco et al. [85], the cumulative rain fallen before the fire event influences soil moisture and its characteristics (e.g., soil infiltration capacity, land use), which control runoff and sediment concentration.
The CN method used in our approach to estimate surface runoff was selected because of its dependence on soil permeability, land use and AMC [64]. To evaluate the highest potential runoff, a CN value corresponding to the worst hydrological conditions (AMCIII) was assigned to each land use. The combination of the CN-AMCIII values with PMP depth, calculated for a 1-day duration and 1.5-year return period, as assumed by Fox et al. [40], allowed to estimate runoff (see Equations (6) and (7)). The variations obtained in the post-fire runoff between 10% and 50% are consistent with the results of other studies [7,42,86,[92][93][94]. Therefore, this approach considers many factors that determine rainfall erosivity. Nevertheless, improvements are plausible by considering other rain durations or by calculating intensityduration-frequency curves.

Erosion
In the erosion module, the erosion values (t/ha) following a PMP for a 1-day duration and 1.5-year return period were calculated in a pre-fire condition and in the two post-fire scenarios-NDVI and RdNBR. Many studies have modeled the erosion assessment, but only 1.4% include an assessment of wildfire effects [95]. This explains the difficulty encountered in our research with regard to sediment concentration data in post-fire scenarios. Furthermore, sediment concentration is site-dependent [42,89,[96][97][98]. For these reasons, data relating to pre-fire conditions but referring to green areas in Apulia [78] were used in our combining approach. This represents a limitation of the model, which may be overcome through field studies following fire events.
The use of the SWAT model sediment concentration spatialized output subdivided into different slope classes, as previously described, allowed to consider slope, which appears to be one of the factors that most influence erosion [99]. In fact, slope is one of the main factors in empirical formulas used for the estimation of erosion [37,100].
To compensate the lack of post-fire sediment concentration data, we introduced two factors that take fire severity into account: NDVI SF and RdNBR SF . These two factors derive from two different spectral indices; the former from the NDVI, directly related to the amount of green biomass, and the latter from the RdNBR, which is sensitive to changes in the amount of live vegetation, moisture content, and some post-fire soil conditions without the pre-fire variability in vegetation composition of NBR [56,62]. With these two severity factors, two post-fire scenarios were developed. Despite the different characteristics of the two factors, the results of the two scenarios obtained are comparable, with differences in the order of a hundredth.
The study results show that rangelands and rangelands with bush are the most vulnerable land uses (i.e., greater increase in erosion values) and are therefore exposed to greater risk. These classes, which exhibited fire-prone behavior with high areal percentage under moderate-and high-fire severity conditions regardless of the area occupied in the site (Figure 3), easily remain bare after the occurrence of a fire. This leads to negatively influence soil properties, such as organic matter content and aggregate stability [101]. Hence, the raindrop effect and, consequently, soil erosion increase [102].

General Remarks and Future Work
The model presented in this work was developed to provide a useful tool to easily assess post-fire erosion and prevent potential risk in subsequent rain events. The three modules of the model are structured so they can be implemented in a GIS environment by users with soft skills in hydrology and satellite image elaborations. Moreover, it does not require much input data, unlike the more complex hydrological models. In addition, its framework allows the model to adapt to different environments and different spatial and temporal scales.
Future studies are needed to compare our combining approach with other validated models (such as models that evaluate post-fire erosion using RUSLE, SWAT or WEPP) to assess its level of correctness. After this phase, more studies are needed to create a comprehensive tool which could be more expeditious for the assessment of vulnerable sites once the wildfire season is over and for the immediate planning of burned area restoration. For example, a single factor could be introduced that takes into account the severity of the fire and other characteristics of the site without considering post-fire sediment concentrations. Additionally, other remote sensing indices can be used to create new scenarios improving the tool according to biophysical and environmental conditions of study sites [84,103].

Conclusions
This study has attempted to provide a combination of methods to estimate post-fire soil erosion for the immediate planning of burned area restoration. Our findings in Apulia are of secondary importance to the approach we have devised because the focal point of the present study is to create the basis for a future tool to be adopted by civil protection or landscape managers. Simplicity and accessibility to a wide range of users make this approach of great practical significance. Furthermore, the production of increasingly accurate and precise remote sensing data will only further improve the use of the tool in supporting policy and management decisions for each fire season at different scales. Once the tool identifies burn sites with the highest post-fire soil erosion, it can be helpful in terms of landscape management for two important reasons. The first deals with the prioritization of soil bioengineering interventions for reducing soil erosion as well as for bolstering protective actions, and the second concerns the budgeting of human and financial resources.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12081105/s1, Table S1: Dates of Sentinel-2 Level-1C products used in the study, Table S2: Areal percentage of land use per study site, Table S3: CN-AMCIII per land use and soil group used for each study site, Figure S1