Multi-Temporal Assessment of Soil Erosion After a Wildfire in Tuscany (Central Italy) Using Google Earth Engine

Abstract

The Massarosa wildfire, which occurred in July 2022 in Northwestern Tuscany (Italy), burned over 800 hectares, leading to significant environmental and geomorphological issues, including an increase in soil erosion rates. This study applied the Revised Universal Soil Loss Equation (RUSLE) model to estimate soil erosion rates with a multi-temporal approach, investigating three main scenarios: before, immediately after, and one-year post-fire. All the analyses were carried out using the Google Earth Engine (GEE) platform with free-access geospatial data and satellite images in order to exploit the cloud computing potentialities. The results indicate a differentiated impact of the fire across the study area, whereby the central parts suffered the highest damages, both in terms of fire-related RUSLE factors and soil loss rates. A sharp increase in erosion rates immediately after the fire was detected, with an increase in maximum soil loss rate from 0.11 ton × ha⁻¹ × yr⁻¹ to 1.29 ton × ha⁻¹ × yr⁻¹, exceeding the precautionary threshold for sustainable soil erosion. In contrast, in the mid-term analysis, the maximum soil loss rate decreased to 0.74 ton × ha⁻¹ × yr⁻¹, although the behavior of the fire-related factors caused an increase in soil erosion variability. The results suggest the need to plan mitigation strategies towards reducing soil erodibility, directly and indirectly, with a continuous monitoring of erosion rates and the application of machine learning algorithms to thoroughly understand the relationships between variables.

1. Introduction

Wildfires represent one of the main causes of environmental deterioration since they affect terrestrial ecosystems and tend to increase geomorphological degradation processes [1], especially in the Mediterranean region [2,3,4]. The current context of climate change (CC) has clear implications in wildfire occurrence in the Mediterranean area for several reasons [5]. The frequency and intensity of heatwaves will tend to increase, leading to a greater likelihood of extreme wildfire events [6]. Regarding the vegetation, a higher number of warmer and drier summers results in greater fuel availability and low moisture content, resulting in an increased fuel flammability and thereby increasing the number of fires [7]. Annual drought periods are strictly correlated to fire occurrence since a direct correlation was found, for instance, in Greece, where two 30-year periods were compared within the second half of the 1900s [8]. Additionally, the CC-related rainfall events with anomalous patterns can result in a lethal mix that boosts the impact of soil erosion phenomena [9]. Fires not only affect the vegetation cover but also alter the physical and chemical properties of soil [10,11], contributing to an increase in its susceptibility to degradation processes like erosion [12]. Although soil erosion is a natural process involving the mechanical removal of soil particles from slopes, primarily due to water runoff, it can become hazardous if the rate of soil loss exceeds the natural soil formation rate [13] or when it triggers other types of geomorphological processes (e.g., landslides). After a wildfire, a significant amplification of erosion intensity is observed since the combustion of vegetation cover leaves the soil without any protection from the actions of water and wind [14]. Several studies have demonstrated the increase in post-fire erosion rates compared to pre-fire conditions, also linked to a reduction in the soil infiltration capacity due to the formation of a hydrophobic layer [15], combustion of organic matter [1], deterioration of soil structure and modification of porosity [16]. The occurrence of heavy rainfall immediately after a fire worsens the general condition of the soil due to nutrient depletion, hindering plant regrowth [17]. Understanding the relationships between wildfires and soil erosion, even with a multi-temporal perspective, is crucial for developing effective land management aimed at mitigating the long-term impact of these events, in particular, in Mediterranean areas such as Italy, where many mountainous and hilly areas are present and complex rainfall patterns are often found, leading to soil erosion processes and landslide events.

Empirical methods are the most used approaches to assess soil erosion rates thanks to their simplicity and flexibility, through the application of the Universal Soil Loss Equation [18] and the Revised Universal Soil Loss Equation [19]. These approaches are applied—particularly nowadays with the improvements of geospatial tools—at the basin scale [20,21], the regional scale [22], and the continental scale [23]. Their ease of application is increased if applied with freely accessible data, allowing them to rapidly assess soil erosion rates over areas of different sizes. Google Earth Engine (GEE) [24] represents one of the most famous platforms for cloud-based geospatial analysis in soil erosion assessment, characterized by a petabyte database of geospatial datasets. Indeed, since its release in 2010, the publications related to soil erosion assessment using the GEE platform have nearly constantly increased from 2016 to 2024. The reduced computational speed also allows for a multi-temporal modeling of environmental parameter changes during seasons or years [25,26], or after forcing events (e.g., floods [27], wildfires [28], landslides [29]).

Although the use of cloud computing is widespread in Italy for environmental tasks (e.g., for gully erosion [30], for conservation agriculture [31], and for archaeological preservation from extreme events [32]) and also for predicting and mapping potential fire severity for risk analysis [33], no particular studies are focused on multi-temporal, post-fire soil erosion assessments in Italy or even consider a mid-term perspective of erosion risk using Google Earth Engine. Therefore, this work aims to overcome this gap by focusing on the multi-temporal assessment of soil erosion rates before and after a wildfire occurs in a small municipality located in Northwestern Tuscany (Italy) using the potentiality of Google Earth Engine in terms of the ease of analyzed time span variation (to assess the soil erosion risk at different time scales) and its complete availability of freely accessible petabytes of data. The Massarosa wildfire is considered as one of the major recent forest fires that has occurred in Italy in recent years. This fire caused damage due to its large spreading area, risking several elements, from the mountainous area down to the coastline. Due to its severity, three scenarios were analyzed by applying the RUSLE approach in order to assess the variation in soil erosion rates after different time spans and to analyze the contribution of vegetation recovery and soil erosion propensity in soil erosion rate changes: (i) before the wildfire, (ii) immediately after the event, and (iii) one year after the event. The RUSLE approach was chosen as opposed to others because of its simplicity, ease of application, and reproducibility that allow for a direct comparison between the three analyzed scenarios and other studies across the Mediterranean area.

2. Study Area

2.1. General Setting

Massarosa is an Italian municipality located in the Northwestern Tuscany region (Lucca province) and it is characterized by different landscapes such as plains, wetlands, and hills, situated in a transition zone between the coastline (Versilia) in the west, the Apennine reliefs, and the Apuan Alps mountain chain in the east and north, respectively, at an altitude that varies from 0 m to 450 m above the sea level. Massarosa presents a typical Mediterranean climate, characterized by hot, dry summers and mild, wet winters with a mean annual precipitation ranging from 870 to 900 mm/year.

From a geological perspective, the hilly slopes are constituted by sandstone bedrocks (Macigno Formation, MAC in Figure 1b) formed during the Late Oligocene to Early Miocene era as a result of turbidite sedimentation, characteristic of deep marine environments during the Appennine belt formation [34], laying over terrigenous formations (Scaglia Toscana Formation, SCA in Figure 1b) and pelagic calcareous rocks (Serie Toscana group, ST in Figure 1b). The land cover (from Corine Land Cover, 2018 update) is mainly composed of forest areas (i.e., the Mediterranean woods and bush forest) in the hilly parts, and agricultural and agroforestry areas as well as urban areas and wetlands in the flat zones, as shown in Figure 1c. Regarding the soil type distribution, the hilly areas of Massarosa are mainly characterized by slightly deep soil associations with an A-Bw-C horizons sequence characterized by a well-drained structure with predominately loamy sand and loam textures (PEL1 facies), less deep soils associations with an A-Bw-R horizons characterized by a well-drained structure with gravelly and pebbly loam textures (SFC1 facies) and slightly deep soil associations with an Ap-C-Cr horizon sequence, well drained with a loamy texture (GCC1 facies) [35]. Near the urbanized areas, deep soils with an Ap-Bw-Bg horizon sequence that are well drained and with loamy and sandy textures (FAB1 facies) are present (Figure 1d).

2.2. The 2022 Wildfire

The Massarosa wildfire was a forest fire that occurred in the mountainous area near the main town center that started on the afternoon of 18 July 2022 and ended on 22 July 2022. Due to the extent of the forest fire, ashes rapidly reached the seaside (approximately 10 Km away from the fire area) and the local “A11” highway section was closed. After two days, a fire-front of several kilometers long was still active with more than 60 hectares burned. Accordingly, the Copernicus Emergency Management Service (CEMS) was activated on 20 July 2022 by the Italian Civil Protection to map the extent of the wildfire, assess the most affected areas, and evaluate the damage. At the end of the wildfire, 836.7 hectares of land was burnt [36,37].

3. Materials and Methods

3.1. Google Earth Engine and Data Sources

Google Earth Engine is a cloud-based platform consisting of a multi-petabyte analysis-ready data catalog and high-performance computational service. The data catalog hosts a large amount of free available geospatial datasets including Earth Observation images from different satellites (e.g., Landsat, Sentinel constellations) and weather, climate, land cover, topographic, and many other types of datasets. The GEE platform can compute geospatial analysis at the planetary scale thanks to the large parallel processing system that subdivides computations, whereby users can gain access to the data and run scripts with the JavaScript programming language using a web-based application programming interface (API) and associated interactive development environment (IDE) [24]. For this study, datasets from different sources within the GEE catalog were used.

As multispectral satellite datasets, three Sentinel-2 L2A (atmospherically corrected) images [38] were collected by filtering the GEE collection on the period of interest and setting a suitable cloud cover percentage threshold (30% for this study). In particular, an image before the wildfire occurrence (acquired on 15 February 2022), an image immediately after the end of the wildfire (acquired on 27 July 2022), and an image after one year from the wildfire occurrence (acquired on 15 July 2023) were used. The topographic variables were extracted from the NASA Shuttle Radar Topography Mission (SRTM) dataset [39]. This dataset provides a digital elevation model on a global scale at a resolution of 1 arc-second (approximately 30 m). To obtain rainfall data, the Climate Hazard Center for Infrared Precipitation (CHIRPS) dataset was used [40] (0.05 degree spatial resolution, approximately 5 m) while soil data (clay and sand percentage, organic matter content, and texture classes) were retrieved, with a 250 m spatial resolution, from the OpenLandMap soil property database available on GEE [41].

3.2. Burned Area Detection and Burn Severity Assessment

The burn area detection and burn severity assessment represents a mandatory propaedeutic analysis in wildfire case studies. This analysis can be performed through the use of the Normalized Burn Ratio (NBR),which is a normalized difference index that can be calculated from a multispectral image with the following Formula (1):

$$\mathit{NBR}={\displaystyle \frac{(NIR+SWIR)}{(NIR-SWIR)}}$$

(1)

where NIR represents the Near-Infrared spectral region and SWIR represents the Shortwave Infrared spectral region. This index is commonly used in burned area detection studies by analyzing the difference in reflectance of vegetated areas (higher reflectance in NIR spectral region than SWIR spectral region) and fire-affected areas (higher reflectance in SWIR spectral region than NIR region). Consequently, high NBR values in (1) represent healthy vegetation, while recently burnt areas show low NBR values. To assess the burnt area extent, the difference between pre- and post-wildfire NBR images ($\Delta NBR$) was calculated following (2) [42,43]:

$$\Delta NBR={\mathit{NBR}}_{prefire}-{\mathit{NBR}}_{postfire}$$

(2)

The higher the values of $\Delta NBR$, the more severe the fire damage. The United States Geological Survey (USGS) proposed a table (

Table 1. Burn severity classes according to $\Delta NBR$ values (from USGS).

Severity Level	$\mathsf{\Delta }\mathit{NBR}$ Values
Enhanced regrowth	<−0.101
Unburned	−0.100 to 0.99
Low severity	0.100 to 0.269
Moderate severity	0.270 to to 0.659
High severity	>0.660

) in order to help in classifying different levels of burn severity.

In this work, the Sentinel-2 bands number 8 (NIR) and number 12 (SWIR) were used to calculate the NBR index for the three scenarios, and $\Delta \mathrm{NBR}$ was evaluated in two situations: after five days from the wildfire occurrence (difference between 15 July 2022 and 27 July 2022 acquisitions) and after a year from the wildfire occurrence (difference between 15 July 2022 and 15 July 2023 acquisitions) and was reclassified using Table 1 to assess the burn severity for the two cases. The actual burned area was detected by masking out all the pixels with $\Delta NBR$ below 0.100.

3.3. RUSLE Parameter Estimation

In this study, the RUSLE model [18,19] was applied to predict the soil erosion rates in the three scenarios by considering the change in the parameters in response to the wildfire. The RUSLE model assesses the soil erosion rates through (3):

$$A=R\ast K\ast LS\ast C\ast P$$

(3)

where A is the soil erosion rate [ton × ha⁻¹ × yr⁻¹], R is the rainfall erosivity factor [MJ × mm× ha⁻¹ × h⁻¹ × yr⁻¹], K is the soil erodibility factor [ton × h × MJ⁻¹ × mm⁻¹], LS is the slope and length factor [dimensionless], P is the support practice factor [dimensionless], and C is the cover factor [dimensionless]. All the RUSLE parameters were evaluated on a pixel-based approach.

3.3.1. Rainfall Erosivity Factor (R)

The rainfall erosivity factor represents the erosive energy of the rainfall. In this study, it was calculated following (4) (modified from [44]) using the CHIRPS dataset:

$$R={\displaystyle \frac{(1163+4.9\ast H-35.2\ast NRD-0.58\ast q)}{100}}$$

(4)

where H [mm/yr] is the mean annual precipitation, q is the elevation at the pixel level (using SRTM DTM), and NRD (number of rainy days) is the mean number of rainy days per year. In this work, a rainy day is considered as a day with more than 1 mm of rainfall. For the first (pre-fire condition) and second (immediately after-fire condition) scenarios, the previous 10 years (2012–2022) were taken into account to calculate the mean annual precipitation. For the third scenario (one-year after-fire condition), the total annual precipitation (from July 2022 to July 2023) was considered.

3.3.2. Soil Erodibility Factor (K)

The soil erodibility factor represents the propensity of soil particles to be removed and transported by water runoff in terms of structure, permeability, and organic matter content. This parameter is strongly influenced by wildfires since the fire modifies these properties, leading to an increase in the K factor values. The K factor is calculated by (5) [18,19,45] using data from OpenLandMap [41]:

$$K={\displaystyle \frac{2.1\ast {10}^{-4}\ast M\ast (12-OM)+(3.25\ast (s-2))+(0.250\ast (p-3))}{100}}\ast 0.1317$$

(5)

where M is a particle size parameter (dependent on the percentage of sand and clay), OM is the organic matter content, s and p are the soil structure code and permeability class, respectively (evaluated from the soil textural classes from OpenLandMap dataset) (

Table 2. Soil textural classes and corresponding soil structure codes.

Soil Texture Class	Soil Structure Code
Sa, LoSa, SaLo	1 (very fine granular)
SaCl, SaClLo, Lo, SiLO, Si	2 (fine granular)
ClLo, SiClLo	3 (medium or coarse granular)
Cl, SiCl	4 (blocky, platy, or massive)

and

Table 3. Soil textural classes and corresponding soil permeability classes.

Soil Texture Class	Soil Permeability Class
Sa	1 (fast and very fast)
LoSa, SaLo	2 (moderately fast)
Lo, SiLo, Si	3 (moderate)
SaClLo, SaCl	4 (moderately slow)
SiClLo, SaCl	5 (slow)
SiCl, Cl	6 (very slow)

) [46].

In the tables above, the abbreviations and their full forms are as follows: Sa: sand; LoSa: loam–sand; Si: silt; SaLo: sand–loam; SiLo: silt–loam; Lo: loam; SaClLo: sand–clay–loam; SiClLo: silt–clay–loam; ClLo: clay–loam; Cl: clay; SiCl: silt–clay; and SaCl: sand–clay.

In this study, to assess the soil erodibility values after the wildfire, K factor original values were multiplied by factors of 1.6, 1.8, and 2 based on the burn severity estimated from $\Delta NBR$ values (low, medium, and high, respectively) (

Table 4. K factor multiplication parameters dependent on burn severity level.

Severity Level	Multiplication Parameter
Low burn severity	1.6
Medium burn severity	1.8
High burn severity	2.0

), following Terranova et al. [44] and Lanorte et al. [47].

3.3.3. Length and Steepness Factor (LS)

The LS factor represents the morphometric and hydorlogic-related parameters that are involved in soil erosion phenomena (e.g., slope, aspect, flow accumulation). For this study, the procedure proposed by Elnashar et al. [46] was followed (using Equations (6) and (8)–(10)) and

Table 5. Steepness factor (S) values based on slope ranges.

Slope Ranges	S values
<5°	S = 10.80 × $sin\left(slope\right)+0.03$
$5^{\circ }$–${10}^{\circ }$	S = 16.80 × $sin\left(slope\right)+0.50$
>10°	S = 21.91 × $sin\left(slope\right)+0.96$

from Wischmeier et al. [18], Desmet et al. [48], and Renard et al. [19].

$$F=|sin(a\left)\right|+|cos(a\left)\right|$$

(6)

$$F={\displaystyle \frac{\frac{sin\left(s\right)}{0.0896}}{3\ast (sin{\left(s\right)}^{0.8}+0.56)}}$$

(7)

$$m={\displaystyle \frac{F}{F+1}}$$

(8)

$$L={\displaystyle \frac{{(FA+scal{e}^2)}^{m+1}-F{A}^{m+1}}{scal{e}^{m+2}\ast X^m\ast 22.{13}^m}}$$

(9)

$$LS=L\ast S$$

(10)

where a is the aspect direction, s is the slope in radians, F is an intermediate factor, FA is the flow accumulation (from Lehner et al. [49]), and m is the slope length exponent (22.13 is the standard slope length exponent).

3.3.4. Cover Management Factor (C)

The cover management factor represents the surface cover and cover management contribution to soil erosion processes, since a higher vegetation cover and appropriate crop management reduce runoff and soil loss [50]. To rapidly estimate the C factor, spectral indices like NDVI (Normalized Difference Vegetation Index) can be used [51,52]. The NDVI index was calculated as the normalized difference between the NIR band (band 8) and red band (band 4) in Sentinel-2 images in the three scenarios considered. For this work, to evaluate the C factor from the NDVI, an equation derived from Van der Knijff et al. [13], Tamene et al. [53], and Elnashar et al. [46] was used (11):

$$C=exp(-2.5\ast {\displaystyle \frac{NDVI}{(1-NDVI)}})$$

(11)

3.3.5. Support Practice Factor (P)

The P factor represents the contribution of agricultural conservation practices to reduce runoff and erosion (e.g., terracing, contouring) [18,19]. In this work, since no information on specific supporting practices were available, a constant value of 1 was assigned to the P factor for the entire area to consider the worst scenario.

4. Results

In Figure 2, the burn severity maps calculated from the $\Delta \mathrm{NBR}$ values are depicted. Figure 2a represents the “five days after the fire” scenario while Figure 2b represents the “one year after” scenario. From Figure 2a, the impact of the wildfire is clearly visible since most of the study area was classified with either a medium or high burn severity, while only a few areas were considered to have been slightly affected by the fire (low burn severity) and just a few spots remained unburned. After a year (Figure 2b), there was a shift towards a less serious scenario, since a decrease in highly damaged areas is evident. However, some sectors classified as high-burn-severity-affected are still present within the area.

In Figure 3, the R factor maps are presented, while the LS factor map is presented in Figure 4. Both of the R factor maps are characterized by a small value range, with maximum values located in the southern parts of the study area and lower values in the central parts. In the one-year-after scenario, an increase in the absolute value range is clearly visible.

The LS factor presents very low values with the maximum value along the higher and steeper parts and the minimum along flat areas.

The factors that were affected by the wildfire effects (K and C factors) are presented in Figure 5 and Figure 6. Considering the pre-fire situation (Figure 5a), the original values of the K factor present low values (from 0 to 0.024 ton × h × MJ⁻¹ × mm⁻¹) and a narrow range between the minimum and the maximum. The higher values are constrained in the lower parts, while the northern part presents lower values of the K factor. Immediately after the wildfire (Figure 5b), a sharp increase in the K factor values is easily noticeable, as it is two times higher than the baseline (upper limit of 0.043 ton × h × MJ⁻¹ × mm⁻¹). The effects of fire are clearly visible in the central sector, near the western boundary of the area, where an increase from 0.02 ton × h × MJ⁻¹ × mm⁻¹ to 0.04 ton × h × MJ⁻¹ × mm⁻¹ was detected. However, the spatial distribution across the study area remains nearly the same as the before-fire condition. One year after the wildfire (Figure 5c), a slight decrease in the K factor values is visible with a reduction in the upper limit of up to 0.039 ton × h × MJ⁻¹ × mm⁻¹. In this scenario, a change in the spatial distribution is detected with a general increase in low K factor values across the area. The boxplots (Figure 5d) highlight a clear difference in K value distributions in the three scenarios. Before the wildfire, there was a relatively low median K factor value with a narrow distribution and a few outliers below the lower bound. The immediate effects of the wildfire show a sharp increase in the median value and a wider distribution, indicating an increase in the K factor variability across the study area as well. After one year, the effects of the wildfire are still present. A high median value of K is visible even though it is slightly lower than the previous scenario, with a wide distribution indicating a variability in soil recovery.

In Figure 6, the variations in the C factor in the three scenarios are presented. The C factor values before the wildfire (Figure 6a) shows a nearly constant spatial distribution of zero values with a few scattered hotspots. Five days after the fire occurrence, the Figure 6b shows a sharp increase in the values of the C factor with a maximum value more than three times higher than the the baseline (from 0.24 in the baseline to 0.79 after the event). A general increase in the C factor spatial distribution is evident, with several peaks of high C values in the central (0.77) and western parts (0.67) of the study area. The mid-term C factor map (Figure 6c) shows a comeback to nearly pre-fire condition for most of the extent of the area that presents values near 0. The maximum values of the C factor (0.34) are constrained into sparse spots that are, in particular, located in the central part, near the western boundary. As can also be seen in Figure 6d, before the wildfire, the C factor distribution presents values clustered to zero with minimal variability and some outliers above the upper bound, indicating a good and homogeneous vegetation cover. Five days after the fire, a clear increase in the C factor value is visible, both regarding the median value and the variability, indicating a reduction in the vegetation cover. Moreover, the high variability and the presence of very-high-value outliers suggest a differential vegetation response to the fire’s impact. In the long term, the C factor returns to general low values (with a median close to zero and a reduced variability but still higher than pre-fire condition), meaning that most of the vegetation has been recovered. However, the presence of numerous outliers confirms the presence of large spots with still sparse and unhealthy vegetation. In Figure 7, the soil loss maps for the three scenarios are presented.

Figure 7 shows the multi-temporal assessment of soil erosion rates using the RUSLE approach in the three scenarios considered. The baseline situation, before the wildfire occurrence, represents a normal condition with very low soil erosion rates (overall mean value at 7.16 × 10⁻⁴ ton × ha⁻¹ × yr⁻¹ and maximum at 0.11 ton × ha⁻¹ × yr⁻¹) and a constant near-zero value distribution in most of the study area. Immediately after the wildfire, a rapid intensification in the soil loss values is observed. The overall mean value increased to 0.19 ton × ha⁻¹ × yr⁻¹ and the maximum value of soil loss reached 1.29 ton × ha⁻¹ × yr⁻¹, indicating an increase in mean values of 264% and more than 1000%, considering the maximum values. A clear change in the spatial distribution of the soil erosion rates is visible, with maxima detected in the central (1.29 ton × ha⁻¹ × yr⁻¹), south (1.08 ton × ha⁻¹ × yr⁻¹), and western (1.23 ton × ha⁻¹ × yr⁻¹) parts of the study area (in accordance with the changing of the K and C factors). After one year, the overall mean soil erosion rate decreases to 0.02 ton × ha⁻¹ × yr⁻¹ (maximum value equal to 0.74 ton × ha⁻¹ × yr⁻¹) with a decrease in mean values of 89% and 42%, considering the maximum values with respect to the immediate post-fire scenario. However, high soil-loss-rate values are still detected in the central part of the study area that is considered the most problematic zone, while values near zero are detected in the rest of the study region (and also in most of the western and southern parts of the area previously affected by the sharp increase in soil loss rates). Nevertheless, considering the mid-term scenario, both the mean and the maximum soil loss values after one year still have higher differences than the baseline situation (more than 26% higher and more than 570% higher, respectively). The analysis of the boxplot distributions (Figure 7d) highlights a low median value and a low variability of the erosion rates before the event, with some outliers presenting high values. Immediately after the fire, despite a sharp increase in median values and a high number of outliers with low erosion rates, no signs of high variability were detected, indicating a low non-uniformity of the fire’s impact, also in terms of erosion rates across the study area. In contrast to mid-term distributions of C and K factors, the soil erosion rates after one year still have higher median values than the baseline conditions, and the main feature observable is an increase in the rates’ variability with respect to the previous scenarios, thereby suggesting a strong differential spatial recovery of the soil loss rates within the study area.

To better understand the change in soil erosion rate distribution, a Kernel density estimation for the three scenarios was evaluated (Figure 8). The use of a logarithm of soil erosion rates allows us to highlight the differences between the three scenarios. Upon analyzing the graph in Figure 8, the differences between pre-fire and immediately after scenarios are focused on a general increase in the absolute values and an increase in the area-under-the-curve width after five days (density peak from 0.7 to 0.25), indicating a slight change in data variability. The mid-term scenario (one year after) presents highly dispersed data compared to the previous scenarios that virtually cover the whole value’s range between the baseline and immediately after distribution, indicating a time-dependent response typology.

Figure 9 presents the differences in soil erosion rates evaluated among the three scenarios and allows us to identify the change in soil erosion rates both from a spatial distribution across the study area (maps on the left side) and the distribution shapes (violin plot on the right side). Figure 9a (left) shows the distribution of the differences in terms of soil erosion rates between the immediate impact of the wildfire and the baseline condition. The main hotspot that can be seen is located in the central part of the area near the western boundary (differences in soil erosion of 1.29 ton × ha⁻¹ × yr⁻¹) and a minor hotspot located in the northwestern part. The rest of the study area shows lower values of soil loss differences with minimum values located in the north zones. The Kernel density estimation plot (Figure 9a, right) highlights a wide distribution of the difference with a maximum of around 0.1, meaning that many parts of the area present nearly constant values before and after the event. However, the skewness of the data towards positive values and the long tail confirm the presence of localized maxima across the study area up to 1.2–1.4 ton × ha⁻¹ × yr⁻¹. In Figure 9b (left), a strong decrease in absolute values is detected. Indeed, most of the area presents lower values of different soil erosion rates compared to the five-days-after scenario, indicating an ongoing recovery of the study area. The Kernel density plot (right side) highlights a general decrease in the soil erosion rate differences, centering near 0 ton × ha⁻¹ × yr⁻¹, and data variability, indicating a reduction in the overall magnitude of the event even with a wider distribution of soil erosion rate values across the area (Figure 8).

5. Discussion

By the $\Delta NBR$ calculation, the burned area detected (area with $\Delta NBR$ > 0.1) was 8.89 Km² (889.65 hectares) wide, in accordance with the final value provided from [37].

The wildfire had different impacts across the study area, highlighting zones affected by distinct burn severities. In particular, the central part is detected as the most affected in both post-fire scenarios, resulting in the most critical as well due to higher value distributions of fire-independent factors (R and LS). This differential behavior is obviously reflected in the post-fire assessment of K and C factors that experienced higher values in the central zones with respect to the whole burned area and an increase in the data variability immediately after the wildfire. Focusing on the “one year after” distribution, a shift between the two factors is visible since the K factor still has a higher variability distribution than the C factor. This result suggests a nearly total and homogeneous recovery of the vegetation, lowering the C factor to almost constant values across most of the area. On the other hand, the soil erodibility presents higher variability due to its direct relations with the burn severity and the physical fire-related effects as a result of the hydrophobic layer formation, an increase in sealing phenomena, decrease in porosity, and organic carbon removal [1,16]. It is important to highlight that the post-fire assessment of the K factor is based on an empirical method that was already used by Lanorte et al. [47] in order to maintain the procedure of K factor calculation to be as easy as possible.

The resulted soil erosion rates are essentially low in all the three scenarios, with mean annual values of 7.16 × 10⁻⁴ ton × ha⁻¹ × yr⁻¹, 0.19 ton × ha⁻¹ × yr⁻¹, and 0.02 ton × ha ⁻¹ × yr⁻¹ for the pre-fire, five-days-after, and one-year-after scenarios, respectively. These values seem to be within the same order of magnitude with those in other studies of similar regions, with similar sizes among the Mediterranean area [54,55,56]. Conversely, other studies applied in Mediterranean areas detected substantially higher mean erosion values [57,58]. It is worth noting that these studies were applied at higher scales (catchment scale) than this study, leading to an increase in overall mean soil erosion rates.

When one analyzes a post-fire soil erosion event, it is crucial to consider the “window of disturbance (WoD)” model [59], which describes the time span in which an increase in sediment yield is observed due to the fire’s impact and which tends to decrease and return to background values after a certain period due to fire-related parameter recovery (e.g., vegetation cover). The time span can vary from a few months to 10 years. In the Mediterranean region, the WoD can vary from 2 to 6 years [9,60,61], but it can even reach 25 years [62]. Given the above time periods, a mid-term analysis should be considered within the window of disturbance of the forest fire. However, the decrease in soil erosion rates, especially for the mean values where a decrease of nearly 50% was detected, suggests a faster recovery than other similar study sites along the Mediterranean area (e.g., Kastridis et al. [62]).

From a broader perspective, the soil erosion rates in the study area prior to the wildfire presents mean and maximum values significantly below the long-term threshold for sustainable soil erosion in the Mediterranean (1 ton × ha⁻¹ × yr⁻¹ [63]). Considering the maximum values detected, immediately after the forest fire, the rapid increase in soil erosion rates leads to rates exceeding the threshold, and thus, to an unsustainable soil erosion rate in the central and northern parts. Conversely, after one year, the rates return to below the precautionary threshold, suggesting a sustainable soil erosion rate in the long term. However, if the analysis is focused on the mean values, the sustainable threshold would not be overcome, as it was not detected in either the short-term analysis or the mid-term one. With a deeper analysis, the higher change in data variability after one year indicates a change in the “soil erosion response” after the event, with an initial immediate response reflected in a substantial increase in absolute values (with a spatial distribution related to the fire-dependent factors) and a “longer-term” response that represents a spatially varying feedback dependent on the natural recovery mechanism.

This behavior may be caused by the different assessments of C and K factors in the mid-term analysis and is mainly led by soil erodibility. Their divergence leads to an increase in soil loss variability over long periods as part of the post-fire resilience of the landscape, whereby interacting with vegetation regrowth and soil physical and chemical properties results in localized areas of normal recovery or prolonged degradation. The soil loss variability in the mid-term analysis also introduces a complexity factor in suggesting mitigation strategies to reduce the erosion over time. The main efforts should be focused on the reduction in soil erodibility, both directly or indirectly. For instance, the addition of soil nutrients and organic matter could help in soil property restoration [64,65]. Furthermore, the vegetation recovery status may not be enough to stabilize slopes in post-fire conditions against erosion since some hotspots of higher values are still present over one year. To overcome this limit, native and deep-root vegetation species could be introduced to enhance the anchoring of soil layers [66].

Lastly, mitigation actions could be conducted by also focusing on fire-independent factors, in particular, the LS factor reshaping the steepest slopes (e.g., adding contour terraces). These slope reduction actions can lead to the slowing down of runoff, allowing for an increase in water infiltration and providing an additional mechanical stabilization. The continuous monitoring of soil loss and vegetation recovery [67] should also be the key point for modifying mitigation strategies and adapting them at current erosion rates.

The use of GEE in this study allowed us to simplify the procedures for estimating soil erosion assessments for several reasons, despite the relatively small size of the study area. Primarily, the use of freely accessible data, even if with coarse spatial resolution, enhanced a reproducibility of the study for other areas with larger size or where high-resolution datasets are not available. The multi-temporal analyses were rapidly carried out without the need of a large amount of data storage. Secondly, the use of cloud computing allowed for a dynamic and rapid assessment, even for future perspectives, in monitoring the recovery processes during the WoD period. Furthermore, the RUSLE method represents a consolidated approach to rapidly evaluate soil erosion rates using fixed and well-known formulas calibrated on generalized site characteristics even if they cannot investigate and recognize patterns among the variables and study non-linearity relationships. However, some drawbacks have arisen from this study. First of all, a full validation with field measurements was not carried out even though the comparison with other similar studies [54,55] suggested that a correct estimation of the overall mean erosion rates is significant. Moreover, the coarse resolution of some datasets (e.g., OpenLandMap) may lead to a punctual overestimation or underestimation of the soil loss in certain areas, and the empirical calculation of the K factor introduces a simplification into the K factor post-fire assessment. Finally, the use of the one-year-based R factor for the mid-term analysis may lead to an overestimation of the R values compared to the trend-based R factor.

6. Conclusions

In July 2022, a huge wildfire hit the mountainous area near the village of Massarosa (Northwestern Tuscany, Italy) causing more than 800 to be hectares burned, more than 200 evacuees, and several environmental and geomorphological issues, including an increase in soil erosion phenomena due to the presence of steep slopes, vegetation cover removal, and soil property modification. This study used free-access data, available on the Google Earth Engine platform, to investigate the multi-temporal soil erosion rate assessment before and after the wildfire by considering three main scenarios (before the event, immediately after the fire, and one year after the fire) and using the RUSLE method. The analysis revealed a differential impact of the fire across the different zones within the study area with the central part that exhibits the highest values of fire-dependent factors (K and C) and the subsequent soil loss rates. The background very-low-soil-erosion rates experienced a sharp increase immediately after the wildfire, causing an overcoming of the sustainable erosion threshold of 1 ton × ha⁻¹ × yr⁻¹ in some hotspots. After one year, the soil erosion rates decreased to below the threshold, showing comparable values with other studies in the Mediterranean region. On the other hand, the variability in the distribution increased due to the different mid-term behaviors of fire-related factors. This particular phenomenon suggests that we focus on mitigation strategies that act, in particular, to reduce the soil erodibility both directly and indirectly and to implement a continuous monitoring of soil loss and vegetation recovery that could be crucial even for other geomorphological associated risks. This work represents a solid approach for studying the wildfire erosion impact in Italy using the Google Earth Engine cloud platform due to its petabyte database of freely accessible datasets and rapid multi-temporal assessment of soil erosion rates. Future improvements will be focused on using Google Earth Engine datasets on wider scales and analyzing non-linear relationships more thoroughly between causes and effects by applying machine learning algorithms.