Testing simple approaches to map sediment mobilisation hotspots after wildfires

Keywords: erosion, index of connectivity, LAPSUS, modelling approach, post-fire, RUSLE, sediment connectivity, sediment export.

In terms of disturbance, fire is usually cited as a critical agent of soil erosion and land degradation, responsible for structuring vegetation dynamics in fire-prone areas, like in the Mediterranean basin (Shakesby 2011). Fires can be responsible for degradation of soil structure, an increase in water repellence, and mobilisation of ash and debris, which can have substantial impacts on surface water quality (Basso et al. 2021). In addition, fires can increase connectivity (i.e. the ease with which sediment is transported; Bracken and Croke 2007; Heckmann et al. 2018), not only in burnt areas but also across an entire catchment (López-Vicente et al. 2020, 2021), allowing efficient transmission of fluxes. These potential fire effects on soil and aquatic resources may increase the costs of managing affected water resources (Fernandez et al. 2003), and have created a strong demand for a post-fire sediment loss prediction tool (Vieira et al. 2014; Parente et al. 2022).

Several models are currently being used to quantify fire impacts on erosion and sediment transport (Girona-García et al. 2021). However, they present several limitations for a rapid and detailed assessment of erosion and sediment transport after fires. For example, models such as the Revised Universal Soil Loss Equation (RUSLE) model and the revised Morgan–Morgan–Finney (MMF) model can simulate erosion in high-resolution grids using low computational resources but they only consider soil loss in situ and not sediment transport. Also, both models have an empirical basis, although the conceptual model driving MMF is more complex (Hosseini et al. 2018; Vieira et al. 2018; Parente et al. 2022). From the models that do simulate sediment movement, the LImburg Soil Erosion Model (LISEM) can simulate sediment transport and deposition in a high-resolution grid, and it has a detailed physical conceptualisation of hydrological and erosion processes but requires a large amount of processing time and capacity, and can therefore only simulate individual rainfall events (Wu et al. 2021b; Vieira et al. 2022). The long-term landscape evolution model Landscape Process Modelling at Multi-Dimensions and Scales (LAPSUS) is also physically based, with the added potential to simulate multiple years at the cost of a less-detailed process conceptualisation, but it still requires a large amount of time and processing capacity (Follmi et al. 2022). Other models, such as the Soil and Water Assessment Tool (SWAT) and the Water Erosion Prediction Project (WEPP), are computationally more efficient than LISEM or LAPSUS despite their complex process conceptualisation, but they simulate space based on aggregated elements (homogenous hydrological response units in SWAT, and interlinked hillslopes in WEPP) that limit the possibility to highlight erosion ‘hotspots’ for further intervention (Basso et al. 2021; Fernández and Vega 2018; Nunes et al. 2018). Furthermore, as reported above, the applications of these models require a large amount of data, including topography, land use, soil type, meteorology, and associated model parameters.

An alternative to these data-demanding models is the use of indices of connectivity, which have the potential to express the general properties of the landscape under evaluation and map the potential connectivity between the different parts of a catchment (Borselli et al. 2008). The concept of sediment connectivity has received growing attention in hydrology and geomorphology over recent decades and can be defined as the degree to which a system facilitates the transfer of sediment (Heckmann et al. 2018). From a range of connectivity models and indices, we can highlight the index of connectivity (IC) proposed by Borselli et al. (2008) and the Sediment Export (SE) tool, which can predict erosion and the transport of sediments in a simple way (Sharp et al. 2020). The IC has been applied to assess the spatial sediment connectivity for different topographic and land cover configurations (Cavalli et al. 2013; López-Vicente et al. 2013), including in a post-fire scenario (López-Vicente et al. 2020; Wu et al. 2021a). In addition, different adaptations to the original IC (Heckmann et al. 2018), such as the aggregated index of connectivity (AIC), are being widely used to assess sediment connectivity when applying different land use management and for post-fire scenario analysis (López-Vicente et al. 2020, 2021; González-Romero et al. 2021, 2022). The SE tool is part of the InVEST (Integrated Valuation of Environmental Services and Tradeoffs) Sediment Delivery Ratio (SDR) model, and is being used to assess sediment transport and retention in catchments (Hamel et al. 2015; Griffin et al. 2020). However, to the best of our knowledge, no study has validated this type of tool to see if it works for post-fire scenarios. Therefore, the following question remains unanswered: Do these tools identify post-fire sediment mobilisation hotspots?

Following this line of reasoning, the present study seeks to compare and quantify the relationship between the simulated erosion and deposition resulting from the LAPSUS model with three other different erosion prediction tools, IC, RUSLE model and SE. Here the erosion and deposition simulated by the LAPSUS model, already calibrated for the study area (Águeda catchment) in a recent study (Follmi et al. 2022), is used as the best representation of erosion and sediment transport in that area due to the unfeasibility of collecting field data at such a large scale. In addition, IC, RUSLE model and SE are used because of their simple representation of erosion and ease of application.

The Águeda catchment is located in north-central Portugal, in the Beira Litoral region, and covers approximately 404 km², with elevations ranging from 8 to 1071 m (Fig. 1a). This region has a wet Mediterranean climate, with a wet period in autumn–spring (October–April) and a dry, warm period in summer (June–September) (Tavares Wahren et al. 2016). According to the latest land use information of the Portuguese authorities (Carta de Uso e Ocupação do Solo, COS, of 2018, DGT 2020), forest covers 76% of the Águeda catchment, agriculture covers 10%, and scrub covers 9%. These forests are a mixture of Eucalyptus globulus (67%) and Pinus pinaster (21%), which are common plantation trees in north-central Portuguese forests. As plantation forests, they are highly managed, with parcels of forest terrain nearby and forest roads criss-crossing the hillslopes; road density in forests is comparable to that in agricultural lands (see Supplementary Material S1). This results in a patchy landscape, composed of small forest plots, each differently managed by a small landowner (see Supplementary Material S2). This spatial configuration of the landscape can be visible on aerial photography, e.g. in Google Earth (which only contains high resolution images after 2005) but is not captured by land use maps. For examples of further descriptions of the Águeda catchment landscape and its evolution, please see Boulet et al. (2015), Doerr et al. (1996), Hawtree et al. (2015), Shakesby et al. (1996), and Tavares Wahren et al. (2016). During the period 1980–2017, 206 fires occurred and more than 50 000 ha burnt in total. The spatial distribution of these fire events indicates that they mainly occurred in the northeastern part of the catchment (Fig. 1b). Finally, post-fire erosion in the Águeda catchment has been studied through field measurements (Shakesby et al. 1996; Nunes et al. 2020) and numerical modelling (Ferreira 1997; Nunes et al. 2018; Follmi et al. 2022).

Fig. 1.  (a) Altitude and main rivers of the Águeda catchment; (b) fire occurrence from 1980 to 2017 and main rivers of the Águeda catchment; (c) location of the Águeda catchment within Portugal.

Follmi et al. (2022) used the LAPSUS model to simulate long-term (decadal scale; 41 years) erosion and deposition dynamics after multiple fires in the Águeda catchment using a cell size of 25 m and an annual time step. LAPSUS is a landscape evolution model that has the potential to investigate long-term and large-scale spatial landscape evolution (Schoorl and Veldkamp 2001) by simulating spatially explicit detachment, transport, and deposition of sediments within a regular grid and at the annual scale. Required model input data include topography (digital elevation model, DEM, derived from a 25 m resolution grid), a soil depth map, changing land use, spatial burn severity, and annual time series of rainfall, evapotranspiration, and infiltration (e.g. Baartman et al. (2012a, 2012b)). In addition, the model uses two factor maps representing erodibility (K_es) and sedimentation potential (P_es). To parameterise the model for both unburnt and post-fire years, Follmi et al. (2022) used a combination of the land use map and the burn severity map, with (for each combination of land use and burn severity) values for K_es and P_es based on the USLE-C factor (vegetation cover factor). For details on the parameterisation and calibration of the LAPSUS model, see Follmi et al. (2022). In summary, Follmi et al. (2022) assessed LAPSUS model output using published data for burnt and unburnt conditions inside the Águeda catchment at multiple scales: individual burnt and unburnt plots (16 m²); fields (0.5–6 ha); and one micro catchment (1 km²) – finding an excellent relationship with observed data (R² = 0.94). Additionally, it should be noted that this model application did not account for the landscape patchiness described earlier because the available land use maps do not provide this information.

This study used the LAPSUS model simulated erosion and deposition as a proxy of real erosion and sediment transport patterns (Table 1) due to the lack of spatially explicit field-observed data for multiple years for large parts of the catchment and the unfeasibility of obtaining these data at such a large scale.

Table 1.  Output of each approach.

The IC was first proposed by Borselli et al. (2008). It estimates the potential link between sediment eroded from hillslopes and the stream system derived from landscape information – for example, land use and topographic characteristics (Borselli et al. 2008). This index can be computed as the logarithm between the upslope component and the downslope component of a landscape, which represents the potential for downward routing of the sediment produced upslope and the weighted distance between two points of the sediment route (Borselli et al. 2008). The resulting map (theoretical range from −∞ to +∞) indicates higher and lower connectivity areas. In this study, we quantified the IC using the SedInConnect tool due to its practical use and availability at https://github.com/HydrogeomorphologyTools (accessed 15 January 2022); it also includes additional information and guidelines for using the tool.

In this study, the IC was applied to identify the hotspots of sediment connectivity in different fire scenarios. Here, we used as weight raster the USLE vegetation cover factors (USLE C-factor) according to the land use and burn severity class for each year to assess the relative effectiveness of land cover management systems in terms of soil loss. Also, the cell size was defined considering the 25 m resolution of the DEM. We performed several parallel tests considering several targets, such as (1) the official river network (according to the Water Framework Directive), (2) the main catchment outlet, and (3) the main polygon of the catchment but excluding IC values from valleys as defined by the Weiss classification (Weiss 2001). Because the results of these different targets were inconclusive, we used the main polygon of the catchment as the target in the final simulations. Finally, the IC output (Table 1) is a raster map with the potential connectivity of the Águeda catchment landscape.

The RUSLE model (Panagos et al. 2015b) was used to estimate mean annual soil loss rates (Table 1), considering only soil detachment and not transport or deposition. Using the gridded data of the European Soil Data Centre (ESDAC), the RUSLE model can be computed for each pixel i following the equation:

where R_i is the rainfall erosivity factor (MJ mm (ha h year)⁻¹), K_i is the soil erodibility factor (t ha h (ha MJ mm)⁻¹), LS_i is a slope length–gradient factor (dimensionless), C_i is a cover–management factor (dimensionless) and was calibrated considering Table 2, and P_i is a support practice factor (dimensionless). This model has recently been used to assess post-fire soil erosion (Miller et al. 2003; Fernández et al. 2010; Karamesouti et al. 2016; Vieira et al. 2018; Blake et al. 2020).

Table 2.  USLE C-factor for each land use and burn severity class.

The SE tool (Sharp et al. 2020) is a spatially explicit model that computes net soil loss (Table 1), combining IC and RUSLE model (Eqn 1) for each pixel i and following two steps. Firstly, the model assesses the amount of annual soil loss using the RUSLE model (Eqn 1). Then it computes the sediment delivery ratio (SDR, Eqn 2), which is the proportion of soil loss that reaches the stream and can be computed as follows:

where SDR_max is the maximum theoretical SDR, set to an average value of 0.8 due to the absence of detailed soil information (Vigiak et al. 2012; Hamel et al. 2015), IC is computed as detailed above (section on Sediment connectivity by Index of Connectivity (IC), IC₀ is set to 1 and k is set to 1.5, and they define the shape of the sigmoid function SDR-IC relationship (Vigiak et al. 2012; Hamel et al. 2015). Lastly, SE is a direct function of the soil loss and SDR factor (Eqn 3), assuming that once sediment reaches the stream, it will end up at the catchment outlet (Sharp et al. 2020).

where SE_i is the Sediment Export for each pixel i.

We are aware that this model has some limitations. For example, IC and SDR are functions of the DEM resolution introducing uncertainty in the predictions; however, we expected low sensitivity to the DEM because IC is a ratio (Hamel et al. 2015). Moreover, the RUSLE model does not simulate erosion processes in rocks (López-Vicente and Navas 2009). In addition, due to the presence of anthropogenic infrastructure (e.g. patches and trails), it is possible that sediment transport is affected by that and not just by natural flow paths (López-Vicente and Navas 2009; Hamel et al. 2015). Patches and trails represent 62% of road density in the local forests, so their contribution to soil degradation might not be negligible.

The outputs of the four approaches (Table 1) were compared for 4 years of high fire activity, namely 1985, 2005, 2016 and 2017. In addition, the years 1985 and 2016 represent years of low to moderate fire severity, and the years 2005 and 2017 represent years of moderate to high fire severity. The comparison considered two steps. Firstly, we analysed the spatial variability distribution in terms of similar patterns and location of predicted soil erosion hotspots. Then we computed several statistical measures of association that are detailed below considering a set of percentile’s (P) classes for each output of each approach, namely ≤5P, 5–25P, 25–75P, 75–90P, 90–95P and ≥95P. In addition, we only used LAPSUS model results for net soil erosion, recalculating net deposition values as zero for this purpose; the results of RUSLE and SE cannot be compared with deposition values because the first represents the average annual soil loss and the second represents sediments exported from that specific grid cell, not accounting for sediment inflows from upstream and therefore also not deposition. Finally, IC cannot be directly compared with LAPSUS model results in terms of values (because IC results are dimensionless), but only in terms of hotspot locations.

The statistical measures of association comprise contingency tables, Cohen’s Kappa coefficient (kappa) (Cohen 1960), the Chi-squared statistical test (χ² test) (Kritzer 1978), and Phi coefficient (φ). The contingency tables include the percentage of area that corresponds to each approach percentile class considering the same percentile class of the LAPSUS model. This parameter was computed using the table function in R environment. Kappa measures the level of agreement between the LAPSUS model P class and the P class predicted by the other approaches, and ranges from −1 to 1 (Parente et al. 2019), which can be interpreted according to Table 3. Kappa was computed as follows:

where p₀ is the probability of agreement and p_e is the hypothetical probability of chance agreement between the LAPSUS model class and the same class predicted by the other approaches.

Table 3.  Cohen’s kappa coefficient (kappa) classes interpretation.

The χ² test indicates whether or not a significant association between LAPSUS model P classes and P classes predicted by the other approaches exists. This parameter was computed using the chisq.test function in R environment. In this study, the following null and alternative hypotheses were used:

H0: (null hypothesis) The LAPSUS model class and the class of the other approaches are independent.

H1: (alternative hypothesis) The LAPSUS model class and the class of the other approaches are not independent. (i.e. they are associated).

The results included the χ² test value and a P-value, which indicates if H0 is true or not. In this study, P-value of less than or equal to 0.05 indicates that H0 should be rejected in favour of H1. To confirm χ² test results, φ was also computed. The latter φ is a symmetrical statistic that represents the type of association between two approaches, and it ranges from −1 to 1; these values can be divided into 11 classes (Table 4). φ was computed as follows:

where χ² is the chi-squared statistical test results, and n is total number of pixels.

Table 4.  Phi coefficient (φ) classes.

In this study, the three approaches computed needed several inputs (Table 5), which generally were the same because the main objective was to compare them. The USLE C-factors used (Table 2) were based on: (1) López-Vicente et al. (2021), which comprises the best approximation to the values of Fernández et al. (2010) for high severity fires, and the values of Vieira et al. (2018) for moderate severity fires; and (2) Follmi et al. (2022) for the other cases. In addition, the values of the USLE C-factor were assigned considering the burn severity computed by Follmi et al. (2022) and the burn severity classes proposed by the United States Geological Survey (USGS 2020). Furthermore, the USLE C-factor measures the impact of cropping and management practices on erosion rates (Onori et al. 2006; López-Vicente and Navas 2009). It ranges from 0 (very effective vegetation) to 1 (cover is so low that it does not decrease the erosion). A high-severity burn has a high USLE C-factor due to the loss of canopy cover and surface cover and reduction in surface roughness (Larsen and MacDonald 2007).

Table 5.  Input parameters for the four approaches.

Several annual soil parameters used in the computation of the RUSLE model and SE (Table 4) were extracted from the ESDAC (Panagos et al. 2015b): (1) P-factor is the support conservation practices factor at a regional level (Panagos et al. 2020); (2) LS-factor is the slope Length and Steepness factor (Panagos et al. 2015a); and (3) K-factor is the Soil Erodibility in Europe at a high resolution (Panagos et al. 2012). Also, the computation of the RUSLE model and SE used a single R-factor for the entire study area, which is the rainfall erosivity factor computed for the Campia weather station (altitude = 448 m; latitude = 40.674°N; longitude = −8.217°W; Fig. 1a) by van der Grift (2021); this was done to ensure comparability with the LAPSUS model results, which used a single rainfall value for the entire catchment (Follmi et al. 2022).

This subsection focuses on results for 2016 and 2017, representing 2 years of widespread fire occurrence with moderate and high burn severity (i.e. the 2016 fire mainly had a moderate burn severity – 69% of the burnt area; Fig. 2c, whereas the 2017 fire mainly had a high burn severity – 47% of the burnt area; Fig. 2d), More information about the percentage of burnt area considering each burn severity class can be found on Table 6). Additionally, these 2 years present high values of the USLE C-factor (Fig. 3).

Table 6.  Percentage of burnt area by year considering each burn severity class.

The model outputs for 2016 (Fig. 4a), corresponding to an extensive burnt area in the south-western part of the catchment (Fig. 2c), show that all four tools could locate the same soil erosion hotspots in general terms. However, some points should be mentioned. In the western side of the burnt area, the LAPSUS model does not simulate erosion (green colours in Fig. 4a1), whereas IC shows (red colours in Fig. 4a2) relatively high values there. Also, relatively high IC values (dark red colours) are computed in the stream network. Additionally, this year’s higher values of the IC in the large burnt area (dark red colours, Fig. 4a2) correspond to the high values of the USLE C-factor (dark red colours, Fig. 3c).

Fig. 2.  Burn severity for the four study years: (a) 1985; (b) 2005; (c) 2016; and (d) 2017, based on the dataset of Föllmi et al. (2022). Main rivers of catchment indicated by blue line.

Fig. 3.  USLE C-factor for the four study years: (a) 1985; (b) 2005; (c) 2016; and (d) 2017. Main rivers of catchment indicated by blue line.

Fig. 4.  Outputs of the four toolsapproaches used in this study for the years (a) 2016 and (b) 2017: (a1, b1) LAPSUS (Landscape Process Modelling at Multi-dimensions and Scales) model; (a2, b2) Index of Connectivity (IC); (a3, b3) Revised version of Universal Soil Loss Equation (RUSLE model); and (a4, b4) Sediment Export (SE). All the maps have main rivers of catchment indicated by blue line.

In contrast, the RUSLE model results (Fig. 4a3) show its high dependency on the USLE C-factor (Fig. 3c) over the whole catchment, having high values (red colours) where the USLE C-factor is also high. Additionally, the RUSLE model (Fig. 4a3) produces, in general, similar spatial patterns to those of the LAPSUS model, showing hotspots in valleys (red colours), similar to the IC results. Finally, SE (Fig. 4a4) shows similar spatial patterns to those of IC but with lower values for streams and higher values for valleys. Also, it presents lower values than LAPSUS and RUSLE models (orange to dark green colours).

The results for 2017 (Fig. 4b) also show a general agreement concerning the area of erosion hotspots for the large burnt area at the NE of the catchment (Fig. 2d); there was less agreement for the burnt area of 2016 (Fig. 2c) in the SW that was disturbed, but not yet fully recovered from the last wildfire. However, there are some points worth mentioning. The LAPSUS model (Fig. 3b1) shows high similarity with the spatial distribution of the USLE C-factor (Fig. 3d) both in burnt areas and with the spatial distribution of burn severity (Fig. 2d). IC (Fig. 4b2) results show a similar spatial distribution to those of the LAPSUS model (Fig. 4b1), but present much fewer hotspots (orange to dark red colours) and a stronger signal (dark orange colour) in the stream network. In addition, it presents the hotspots (orange to dark red colours) in the burnt area of 2016 in the higher values (dark red colour) of the USLE C-factor (Fig. 3d) and high severity (Fig. 2d). RUSLE model (Fig. 4b3) results show the same spatial distribution as the burn severity (Fig. 2d) and as the USLE C-factor (Fig. 3d), but with high values in the southwest of the catchment. In addition, the RUSLE model shows lower values (light yellow to dark green colours), whereas the LAPSUS model indicates high values (light yellow to light green colours), especially in the north-east part of the catchment. Additionally, on this same north-east part of the catchment the green dark pattern of RUSLE model suggests a similar pattern of unburnt areas of Fig. 2d. Also, the RUSLE model shows hotspots in valleys (orange to dark red colours) that the LAPSUS model does not show, especially in the region of the 2016 burnt area (orange to dark red colours). Finally, SE (Fig. 4b4) presents patterns similar to the LAPSUS model, but with much lower values (light orange to dark green colours). Additionally, unlike the LAPSUS model, SE does not identify hotspots in the mountains, and SE does not present any similarity with USLE C-factor and burn severity spatial distributions.

The results of the contingency tables indicate a similar relationship between the LAPSUS model and the other tools for all the years (Fig. 5). From these results, we can highlight several points. Only 5% of each tool in the years 1985, 2005 and 2017 coincide with LAPSUS model values ≤5P, and 10% in the year 2016. Additionally, it is noteworthy that some approaches have a better relationship with the LAPSUS model than others, and that this relationship changes with the year. For example, from 50 to 61% of the values of RUSLE model and SE above the 75th percentile for the first 3 years match with those of the LAPSUS model. On the other hand, for the year 2017, 75 and 86% of the values above the 75th percentile of the RUSLE model and SE, respectively, fall in the same places on the map as those of the LAPSUS model. Also, 64 and 66% of the values above the 75th percentile of IC for the years 1985 and 2016, respectively, and 82 and 86% for the years 2005 and 2017, respectively, match with the same class of LAPSUS model.

Fig. 5.  Outputs of the contingency tables between LAPSUS (Landscape Process Modelling at Multi-dimensions and Scales) model percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P) and the Index of Connectivity (IC, left column) percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P), the Revised version of Universal Soil Loss Equation (RUSLE model, middle column) percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P), and the Sediment Export (SE, right column) percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P) for the study years, from upper to lower row line 1985, 2005, 2016, and 2017. These results are shown in stacked area graphs.

When we compare the values above the 90th percentile of all three tools with those of the LAPSUS model, the relationship changes a little. For example, more than 90% of IC values, in the years 1985 and 2005 – and of the RUSLE model in the years 1985 and 2016 – fall into these extreme LAPSUS model values. Additionally, from 88 to 91% of SE values for the first 3 years fall in the same places on the map as those of the LAPSUS model. It should be highlighted that these values are, on average, 30% higher than the ones for the relationship between the approaches and the LAPSUS model for the values above the 75th percentile. Finally, in 2017, 91% of IC values and 83% of RUSLE model values matched the values above the 90th percentile of the LAPSUS model. In addition, 95% of SE values above the 90th percentile match those above the 90th percentile of the LAPSUS model. These values above the 90th percentile of 2017 are from 5 to 10% more than those above the 75th percentile of the LAPSUS model. The difference between 2017 and the other years might be explained due to the presence of vegetation recovering from the fire of 2016.

The statistical measure of association, kappa (Fig. 6), indicates there is fair to substantial agreement (0.21 < kappa < 0.8) between the values above the 95th percentile of all three tools with those of the LAPSUS model. Additionally, for the values from the 5th percentile to the 95th percentile, kappa points to no to fair agreement (0 < kappa < 0.4) between the three tools with those of the LAPSUS model. In addition, the values ≤5P kappa indicates near perfect or acceptable agreement (0.8 < kappa < 1) between the three tools and those of the LAPSUS model.

Fig. 6.  Cohen’s kappa coefficient (kappa) outputs between the LAPSUS (Landscape Process Modelling at Multi-dimensions and Scales) model percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P) and the same percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P) predicted by the other tools (IC, Index of Connectivity; RUSLE, Revised version of Universal Soil Loss Equation model; SE, Sediment Export) for the study years, from upper to lower row 1985, 2005, 2016, and 2017. Horizontal lines indicate thresholds of kappa classes.

Overall, the χ² test (Fig. 7, columns) shows significant differences (P < 0.05 with four degrees of freedom) for all the relationships between the tools and the LAPSUS model, which points to an association between them; however, our χ² test values are much larger than the critical value 9.49 (refer to the table of ‘critical values of Chi-Square’ in Fisher (1992), and is indicated in Fig. 7 as a red line), which means that this association is unlikely to occur. Therefore, the LAPSUS model values differ significantly from the other tools’ values for all percentile classes. Finally, the φ results (Fig. 7, black lines) confirm this result classifying this association as being non-existent or negligible (0 < φ < 0.2).

Fig. 7.  Chi-squared statistical test (χ² test, in columns) and Phi coefficient (φ, in lines) outputs between the LAPSUS (Landscape Process Modelling at Multi-dimensions and Scales) model percentile (P) classes (<5P, 5–25P, 25–75P, 75–90P, 90–95P, and ≥95P) and the same P classes of the other tools (IC, Index of Connectivity; RUSLE, Revised version of Universal Soil Loss Equation; SE, Sediment Export) for the study years, upper to lower row 1985, 2005, 2016, and 2017. Critical Chi-squared value (CV) indicated by red line.

The main outcome of this study is that similar spatial patterns of sediment mobilisation hotspots (in this case, defined as those with mobilisation rates above the 95th percentile) could be observed in all the tools but with somewhat different locations according to the year (orange to dark red colours of Fig. 4, Supplementary Fig. S1). Furthermore, the results of the contingency tables (Fig. 6) suggest that only LAPSUS model values above the 90th percentile are well represented in all the tested tools, which means that all three tools indicate the same erosion hotspots but not the areas of intermediate erosion. This might indicate that they will be helpful when the aim is to locate areas with extreme erosion rates. In addition, this relation varies with fire characteristics (fire location, fire severity, and burnt area) from year to year, which might mean that these tools require calibration according to fire characteristics.

The results offer several clues to the reasons for these discrepancies. An analysis of the contingency tables (Fig. 6) suggests that a large part of the errors is the result of the misrepresentation of areas with LAPSUS model values above the 75th percentile, which can correspond from 30 to 70% of areas with values below the 75th percentile in the other three tools. These results suggest that the distribution of erosion frequencies is much more skewed towards higher values in the RUSLE model, IC and SE than in the LAPSUS model (i.e. erosion hotspots are more concentrated in these simple tools). This could result from the fact that they do not consider sediment transport capacity, thus predicting high erosion values in regions of concentrated flow where runoff transport capacity is usually high. Also, IC and SE predict hotspots along the stream network that are not present in the LAPSUS model. This is a consequence of not considering sediment movement because the stream network indeed has high connectivity. However, the transport capacity of streamflow is low due to shallow local slopes, resulting in high sediment deposition, which the LAPSUS model simulates. In combination, these results suggest that these simple tools considering only detachment work well in predicting hotspots because they are located in areas where erosion is detachment limited. Additionally, this limits their ability to assess areas with intermediate erosion ranges that might be transport limited. The detachment limitation in the upper reaches vs transport limitation in lower reaches of a burnt catchment has also been described by Wu et al. (2021b), who used a Topographic Position Index (TPI) to identify these reaches; the TPI could be combined with the application of IC and SE to identify and remove these areas with potential false positives.

Another interesting result is that there is a substantial agreement between LAPSUS and the other tools in the values above 95P in all years except 2017, and in that year, those differences mainly relate to an area in the south-western part of the catchment that burnt in 2016 and was in a recovery process. This indicates that the ability of these tools to agree with the LAPSUS model is much more significant in the first year after the fire, and that they are less robust when analysing erosion hotspots after some vegetation recovery has taken place. In addition, our results for intermediate erosion rates (values from 5P to 95P) suggest that the tools might not indicate the same type of processes in this range. Also, the near perfect agreement for values below 5P might indicate that the tools agree well where there are hardly any erosion processes. Finally, our analysis of the statistical measures indicates that, overall, the association between LAPSUS and the other tools is unlikely to occur.

Unfortunately, there are no studies comparing these three tools (IC, RUSLE model, and SE) across different burnt landscapes to compare with our results. However, previous studies have analysed these tools’ performance and can provide some indications of why they work. The IC has been tested in burnt areas by Fernández et al. (2020). They identified it as an excellent tool for assessing post-fire sediment delivery using USLE C-factor based on soil burn severity. Although vegetation burn severity was used in this study as a proxy, this was also used in other studies (López-Vicente et al. 2021) and should not be considered a significant limitation. The RUSLE model has been shown to have acceptable performance in small burnt plots by Vieira et al. (2018). In both cases, a good performance for hotspots was already established. To our knowledge, the SE has not been assessed for burnt areas, but because it results from a combination of IC and the RUSLE model, we could also expect that it would provide good results. One crucial point is that all tools used the same approach to calculate the USLE C-factor from burn severity, which probably contributed to the similarity in the results. However, our results point to a more prominent dependance of the spatial distribution of the RUSLE model on the USLE C-factor than of the spatial distributions of the IC and of the SE (Figs 3, 4, Supplementary Fig. S1).

One important limitation of this study is the consideration of the LAPSUS model as a good proxy for real erosion patterns. Being a landscape evolution model, LAPSUS simulates erosion and deposition on an annual timescale and depends, for parameterisation, on (historic) land use and burn severity maps. Such maps may be inaccurate or lack details such as roads or patchy land use patterns. This limitation could potentially be overcome by parameterising vegetation cover using satellite imagery directly (van Eck et al. 2016). However, these relationships are usually not linear (Ayalew et al. 2020), especially in burnt areas (Fernández and Vega 2016), so this might add an unknown amount of uncertainty to the model. Furthermore, the annual time step does not allow the simulation of severe individual events, which could be important erosion agents in burnt areas (Nunes et al. 2020; Wu et al. 2021a). To improve this, other models that include rainfall events might be applied (Wu et al. 2021b), but these models are often more data demanding, as discussed in the introduction.

The performance of the LAPSUS model in simulating the impacts of wildfire on erosion has been assessed for this area at multiple scales, including measured erosion rates at a burnt hillslope (5.9 ha), two agricultural areas (0.63 ha), and turbidity measurements at the outlet of a small headwater catchment of 94 ha (Follmi et al. 2022). In addition, plots (8 × 2 m) were used to calibrate unburnt erosion rates, but the areas of the assessed wildfires were much larger and had a higher variation of burn severity and slope conditions than those for which the model was calibrated. These calibration issues, and (common) lack of calibration and validation data are highlighted multiple times in Follmi et al. (2022) and also apply in this study. Therefore, Follmi et al. (2022) also indicates that the LAPSUS model results are more suited to showing the spatial patterns of erosion and deposition, and give an indication of the order of magnitude of their rates, than to making precise predictions. The use of simulated results is difficult to avoid due to the strong logistical limitations in carrying out erosion assessments for large burnt areas. However, it would be interesting to test these indices for areas where erosion patterns were measured, even if using simplified approaches such as visual surveys. An alternative would be to first test the indices on smaller catchments where field data on erosion and deposition is available, and only then applying them to larger regions. Also, because the patchy local land use pattern was depicted on the land use map for the entire watershed, it was not included in the LAPSUS model; it could therefore be tested in a sub-catchment where data are available, to assess the uncertainty added by this patchiness to the results.

In summary, the results indicate that for years of low to moderate severity (e.g. in the years 1985; Supplementary Fig. S1a, and 2016; Fig. 4a, more than 80% of the burnt area was in low to moderate severity), the RUSLE model (Fig. 4a3, Supplementary Fig. S1a3) gives the best match of spatial erosion patterns with the LAPSUS model (3% more of the values than IC and SE). Conversely, for the year 2005 (Supplementary Fig. S1b), a year of moderate to high severity (86% of burnt area was in moderate to high severity), IC (Supplementary Fig. S1b2) shows spatial patterns that fit with those of the LAPSUS model better (6% more of the values than those of the other tools), and SE (Fig. 4b4) gives the best spatial patterns match with those of LAPSUS model during the year 2017 (4% more of the values than IC; Fig. 4b2, and 12% more of the values than the RUSLE model; Fig. 4b3). However, it must be highlighted that differences between tools are generally around 5%, which indicates that all three tools perform equally well in locating erosion hotspots, defined as regions with erosion above the 95th percentile, in the first year after the fire.

Finally, we want to underline that Follmi et al. (2022) did not assess the correctness of the erosion rates as predicted by the RUSLE model and SE, and that IC does not predict actual rates. However, the hotspots identified in our maps/results can be used to target post-fire intervention actions to control erosion rates. These kinds of actions have had some recent developments; for example: (1) the Portuguese Institute for Conservation of Nature and Forests (ICNF, https://icnf.pt/, accessed 12 February 2022) has been tasked to assess and implement soil conservation measures; and (2) efforts to implement mulching using voluntary work have proved to be a low-cost but effective intervention option (Prats et al. 2022). Also, we believe that these recent developments can be supported using rapid assessment methods such as IC. Finally, this study has focused on assessing how closely the different erosion-prediction tools agree as regards potential erosion hotspots following wildfire. The influence of land management practices lies outside the scope of this study but, at small scales, they could be expected to contribute to erosion (Martins et al. 2013). Whilst the tools outputs would provide the broad picture of erosion hotspots, these small-scale effects might also need to be identified/borne in mind in determining the most effective distribution of resources for mitigating post-fire erosion throughout a large area like that of the Agueda catchment.

This study aimed to identify sediment mobilisation hotspots after fire using three relatively simple erosion prediction tools, namely the RUSLE model, IC and SE, in a large catchment (404 km²). To achieve this, we compared these tools with simulated erosion by the LAPSUS model, a landscape evolution model (already calibrated for the study area) that could simulate annual erosion response in recently burnt areas. As such, two significant results can be highlighted:

1. The three tools tested in this study are suitable for identifying sediment mobilisation hotspots (i.e. with erosion rates above the 95th percentile, in the first year after the fire), but they are not suitable for areas with erosion rates from 5th to 95th percentiles, which could still be necessary.

2. Although the best-performing approach varies according to fire extent and characteristics, the differences between them are minor, and all three perform equally well for hotspot mapping.

The results of this study can be used for post-fire and water contamination risk management because these tools are relatively simple to apply using information that can be quickly obtained after a fire, thus allowing for a rapid prioritisation of areas for emergency post-fire intervention. However, because the LAPSUS erosion modelling results were uncertain and calibrated on the only available sub-catchment observed data, a next step would be to further verify these results with observed erosion and deposition data in other post-fire affected areas.

Supplementary material is available online.

The data that support the findings of this study are available from the corresponding author, Joana Parente, on reasonable request.

The authors declare no conflicts of interest.

This work was financed by the project FRISCO (managing Fire-induced RISks of water quality Contamination (PCIF/MPG/0044/2018)) and with funding attributed to the CE3C research center (UIDB/00329/2020).