Prediction, validation, and uncertainties of a nation-wide post-fire soil erosion risk assessment in Portugal

Wildfires are a recurrent and increasing threat in mainland Portugal, where over 4.5 million hectares of forests and scrublands have burned over the last 38 years. These fire-affected landscapes have suffered an intensification of soil erosion processes, which can negatively affect soil carbon storage, reduce fertility and forest productivity, and can become a source of pollutants. The main objective of the present study is to produce a post-fire soil erosion risk map for the forest and shrubland areas in mainland Portugal and assess its reliability. To this end, the semi-empirical Morgan–Morgan–Finney erosion model was used to assess the potential post-fire soil erosion according to distinct burn severity and climate scenarios, and the accuracy of the predictions was verified by an uncertainty analysis and validated against independent field datasets. The proposed approach successfully allowed mapping post-fire soil erosion in Portugal and identified the areas with higher post-fire erosion risk for past and future climate extremes. The outcomes of this study comprise a set of tools to help forest managers in their decision-making for post-fire emergency stabilization, ensuring the adequate selection of areas for mitigation to minimize the economic and environmental losses caused by fire-enhanced soil erosion.


Materials and methods
In this study, we developed a post-fire soil erosion risk map for Portugal using the revised MMF model, by applying a successful field-based parameterization, conducted in previous studies, that addressed different burn severities and land uses [21][22][23] (Fig. 1). The development of this map required several steps of data processing, to adjust variables from different sources to the same spatial resolution and to the model input requirements. Afterwards, we computed the revised MMF model for the 38 hydrological years between October 1st 1980 and September 30th 2018, using two different rainfall datasets and for three different soil burn severity scenarios. The inputs for the revised MMF model consisted mostly of spatially distributed data characterizing topography, soil properties, and rainfall regime.   Study area. Portugal is located in Southwestern Europe with a total area ca. 900,000 km 2 presenting an extensive history of wildfire recurrence (Fig. 2) 9,27 . The mainland area of Portugal presents an altitude ranging from sea level in the western and southern coastal areas to about 2,000 m a.s.l. in the north-central region. It is located in a transition zone between sub-tropical and mid-latitude climates 8 , where temperate climate prevails with dry and warm summer (Csb) in the north, and with dry and hot summer (Csa) in the south 28 . The type of climate depends on the topographical features of the country and helps explaining the land cover differences among regions (Fig. 3a). Those climate conditions increase the wildfire susceptibility of Portugal, and explains the prominent annual wildfire activity cycle (Fig. 2) 9 .
Spatially distributed datasets. The gridded global atmospheric reanalysis products of the European Centre for Medium-Range Weather Forecasts ERA-Interim and ERA5 datasets were selected to obtain annual rainfall data (Fig. 3). Both datasets were selected because MMF requires as input a continuous series of daily precipitation data 31,32 .
A high-resolution map of slope angle was derived from the Digital Elevation Model that resulted from the Shuttle Radar Topographic Mission (Table 1). The map showed a heterogeneous slope steepness distribution within Portugal, observing the highest values in the north and southernmost regions (Fig. 4a). Soil bulk density  [28][29][30] . Points indicate the location of the validation sites defined in   Table 1.  (Table 1). Land cover national information was obtained from the national inventory of 2018 provided by Direção Geral do Território. It comprises 83 land cover sub-classes but, given the focus on wildfires, we used five (Table 1). Together, these 5 sub-classes correspond to about 3.3 million ha, which comprise 63% of the main forest types of Portugal and all the shrubland area 37 . The forests dominated by eucalyptus (Eucalyptus spp.) and maritime pine (Pinus pinaster Ait.) are mostly located in Central Portugal (Fig. 4b), while the stone pine (Pinus pinea L.) forests are mainly located immediately at the top of south region (Fig. 4b). The shrubland areas are the dominant land cover in the northeastern and southernmost regions of Portugal.
The soil bulk density, soil moisture content, and ERA5 datasets needed to be scaled down to fit the 25 m grid dimension (Table 1). To this end, the Inverse-Distance Weighting method (IDW) was used, which is a spatial interpolation deterministic model that was computed using R gstat package 38 . Reference soil erosion and local precipitation. MMF performance as well as the global rainfall datasets were validated against published field observations. To this end, the SCOPUS dataset was searched to compile all peer-reviewed articles in English in international journals that reported annual post-fire erosion rates in Portugal. The resulting 11 articles included a total of 10 field sites with micro-plot to swale-scale erosion data on 3 to 18 plots per site, giving a total of 20 years of study ( Table 2). The articles also reported the relevant rainfall data for all 10 sites. There was a strong bias in the geographical location of the field sites, the bulk of them being located in north-central Portugal. This bias, however, largely coincides with the strong prevalence of burned areas in north-central Portugal over the past decades (Fig. 2), and also reflects the lack of post-fire erosion studies in the other parts of the country.
Post-fire soil erosion predictions by the MMF model. The MMF model predicts annual runoff and associated soil losses at slope scale 22 , separating the processes in a water phase and a sediment phase 39 (Fig. 1). The first phase simulates overland flow generation as well as detachment of soil particles by rainfall and runoff, while the second phase models the transport capacity of the overland flow to transport the detached soil particles 22,39 .  www.nature.com/scientificreports/ In this study, the MMF model was applied to predict annual soil losses for a 38-year period, considering three scenarios of soil burn severity (low, moderate, and high). The maximum annual erosion rates predicted over this period of time were then used for mapping post-fire erosion risk. The MMF parameters and associated model input datasets were the following: 1. Annual rainfall (mm yr −1 ) and days with rain: derived from the above-mentioned global rainfall datasets; 2. Soil parameters: bulk density (BD, g cm −3 ), soil detachability index (K, g J −1 ), and soil surface cohesion (COH, kPa), were estimated from LUCAS dataset 32 . The effective hydrological depth of soil (EHD, m) was calibrated with field data 14,40 , and previous MMF applications in recently burned areas in north-central Portugal 21-23 ( Fig. 1); 3. Landform: represented as slope steepness (S, °); 4. Land cover parameters: rainfall interception (A) and canopy cover (CC, %), plant height (PH, m) were approximated to 0 given the wildfire effects. The ratio of actual (Et, mm) to potential (E0, mm) evapotranspiration, crop cover management factor (C), and ground cover (GC, %) and estimated from previous MMF applications in recently burned areas in north-central Portugal 21-23 (Fig. 1); 5. Maximum soil moisture at storage capacity: estimated from soil moisture dataset.

MMF uncertainty analysis and validation.
The parametric uncertainty analysis of the MMF erosion predictions was performed with respect to the rainfall input data. To this end, the predictions resulting from the ERA-Interim and ERA5 datasets were compared in terms of total predicted soil losses, magnitude, and spatial discrepancies. In addition, MMF performance was assessed by comparing the predicted and observed post-fire erosion rates for the 10 field sites listed in Table 2. This was done using the rainfall data from the two global rainfall datasets and from the field studies. To quantify MMF performance, the following four indicators were used in line with previous studies 21  www.nature.com/scientificreports/ It should be noted, however, that such metrics are more frequently used within continuous catchment hydrology modelling approaches than for spatially distributed predictions 39 . Despite those limitations, such methodology is still the most suitable one to assess model performance and adaptation.

Results
Post-fire soil erosion risk mapping. The nation-wide application of the MMF model to Portugal revealed substantial differences in terms of predicted soil erosion amounts based on the rainfall dataset choice (Fig. 5). For the entire studied period (1980 to 2018), the highest soil erosion predictions were found for the year 2001, Nevertheless, in both cases the highest post-fire soil erosion risk classes are predominantly found in North-Central Portugal, coinciding with a greater distribution of steep slopes (Fig. 4) and high rainfall amounts (Fig. 3), over areas with shrublands and eucalyptus as dominant vegetation cover (Fig. 4).
When considering the effect of burn severity in our modelling predictions, it is noteworthy that soil burn severity dramatically increases sediment losses, being this aspect more evident for ERA-Interim than for ERA5 (Fig. 5). For low burn severity, predicted soil losses are mostly (99.9%) below the 10 Mg ha −1 yr −1 threshold for ERA-Interim, while for ERA5 this limit decreases to 1 Mg ha −1 yr −1 . For moderate burn severity conditions, the areas with predicted post-fire erosion risk above the threshold of 10 Mg ha −1 yr −1 substantially increased for the North of Portugal using the ERA-Interim dataset, while for ERA5 these rates are still below. For high burn severity, 14% of the target area presents a severe risk of post-fire soil erosion, with rates above 50 Mg ha −1 yr −1 for the ERA-Interim dataset, and predominating 1 to 2 Mg ha −1 yr −1 classes for the ERA5 dataset.

Model uncertainty analysis and validation.
Our results show that, for the high burn severity scenario, soil erosion predictions obtained using the ERA-Interim data set are two orders of magnitude higher than those using the ERA5 data set. Despite the fact that the spatial distribution of rainfall of both datasets presents slight differences (Fig. 4), no major spatial discrepancies were found when predicting post-fire soil erosion (Fig. 6); being the ones resulting from ERA-Interim consistently higher than those from ERA5 because of its almost consistently higher annual rainfall amounts, with the exception of few locations (0.002%).
Overall, the MMF predictions for the field sites presented a poor agreement with the observed average postfire erosion rates, which could be related to their high variability (Table 3; Fig. 7). In terms of all four indicators, MMF performance was not satisfactory for any of the four rainfall datasets according to the thresholds proposed in Moriasi et al. 41 . Nonetheless, MMF performance did differ noticeably between the four data sets. Overall, it was best when using the local annual rainfall data combined with the ERA-interim data (number of raining days), and worst when using the ERA5 dataset. Furthermore, the use of local rainfall data generally improved model performance for the two global rainfall datasets, except in terms of R 2 in the case of the ERA5 dataset. Finally, MMF performance was markedly better for the ERA-Interim than for the ERA5 dataset, in terms of all four indicators. In this way, the ERA5 dataset led to a higher underestimation of observed erosion rates than ERA-Interim, as shown by the negative PBIAS (Table 3) and Fig. 8.

Discussion
This study aims to produce a soil erosion risk map for the wildfire-prone vegetated areas in Portugal 26 by using the revised MMF model. This model has been previously calibrated for Portugal 21,22 and NW Spain 23 and has demonstrated great capacity to simulate annually and seasonally the hydrological and erosive response in recently burned forest areas 21,22,42 . As a final outcome of this study, such post-fire soil erosion risk map for Portugal was possible to be created with the available data, although several important uncertainties were identified during its development. Nevertheless, given the uncertainty analysis and the model performance assessment made with the local validation sites, the high severity scenario obtained using ERA-Interim was selected for the nation-wide post-fire erosion risk map, as its estimates agreed best with the field observations (Fig. 9).
Despite areas at risk of extreme post-fire soil erosion rates (> 10 Mg ha −1 y −1 ) were identified in both north and south of Portugal, the north-central part of the country showed a much higher density of locations with such risk (Fig. 5). This can be justified by the greater eucalypt, pine and shrubland cover areas (Fig. 3), and the highest distribution of steep slopes at these locations. Such result also agrees with previous studies, and despite Panagos et al. 43 used a different rainfall input in the form of rainfall erosivity, their soil erosion risk obtained for Europe by the Revised Universal Loss Equation model (RUSLE) presents a similar spatial pattern to ours 44,45 .
This implies that in order to prevent post-fire soil erosion in recently burned areas, the decision-makers and land managers of the north-central region areas of Portugal need to be aware of such risks. Especially considering that those high-risk areas are recurrently affected by wildfires as the historical data shows (Fig. 2), and that these same locations also provide important ecosystems services to society such as the maintenance of water quality and flood control, which can be severely affected by post-fire soil erosion 16,46,47 .
This post-fire soil erosion risk map (Fig. 9) will help managers in the decision-making process after wildfires, allowing the identification of high-risk areas and thus, of priority for intervention, ensuring the most efficient implementation of post-fire mitigation measures 19 . It should be highlighted however, that the transformation of the predictions obtained in this study to an easily accessible tool for post-fire management would further improve www.nature.com/scientificreports/ decision-making in these high-risk areas, by considering the various scenarios of severity and the specificities of each affected area in a dynamic way. A parametric uncertainty analysis was performed to be used as a guideline in the decision-making process when the purpose is to consider and apply, on time, the proper soil erosion mitigation measures after wildfire. Taking this in mind, the methodology used to create the predictions was practical, scientifically defensible, and robust enough to be applicable. For validation matters, several field data were considered as the best estimation possible of the post-fire soil erosion rates.
The differences between the predictions from ERA-Interim and ERA5 (Fig. 6) underline that different precipitation datasets can contain significantly different spatiotemporal information, which might lead to contrasting soil erosion estimations. Also, the authors acknowledge that these discrepancies might be from the limitations of IDW used to scaled down ERA5, which is in spite of being simple, fast, and easy to compute and to interpret 48 , has some limitations such as: (i) weighting parameters are chosen a priori, and not empirically; (ii) the variances of predicted values in unsampled locations cannot be estimated; and (iii) the IDW standard application assumes that the distance-decay relationship is constant through space. In this respect, it is important to underline the good agreement between total annual mean precipitation results obtained in this study with ERA-Interim for 1980-2018 period (Fig. 4) and the maps published in the Iberian climate Atlas produced with data observed in a large set of weather stations, although for a slightly different (1971-2000) period 49 .
Finally, the authors acknowledge that biases of the reanalysis may affect the agreement between the soil erosion predictions and soil erosion. In addition to that, some limitations were also found concerning the available methodologies to assess model performance, especially for validation purposes. By using Moriasi et al. 41 metrics, only average soil erosion measurements were used, which did not account for the high variability of field results. Notwithstanding, ERA-Interim seems to be the most advisable reanalysis product because it produced the best results (Fig. 7), despite the underestimations found when compared with the field validation data. Which is in line with Belo-Pereira et al. 50 , who have demonstrated that ERA-Interim underestimates precipitation in the rainiest months in Western Iberia and tends to overestimate it in NE and SE regions.
It is therefore advisable that authors consider more than one source of data when creating risk maps to account for additional uncertainties in the predictions; in addition to considering local data when available to better evaluate post-fire soil erosion assessment, as illustrated by the two main sources of uncertainties highlighted in this study (severity and rainfall). Those considerations would eventually translate into more detailed and accurate risk maps and thus, better-informed land management decisions 24 . Identifying all the possible outcomes can help land managers and policy makers to understand that there is always some level of uncertainties associated to post-fire responses, and that those uncertainties should be included in the decision making. In addition to that, it should be highlighted that additional sources of uncertainties other than those pointed out in this study can also be included in such risk assessment. Aspects such as future climate demands, assessment scale, land management operations, or even the inclusion of land uses different to those identified in this study can improve the understanding of the risk at hand; whether the focus is the mitigation of the impacts of the current wildfire, or the prevention of future post-fire impacts, but also, if such reflections are made at local scale or at a political level.
Moreover, it should be highlighted that wildfire occurrence is not exclusive from the Portuguese territory and corresponds to a global environmental problem 51 , and therefore the application of such approach should be considered at the European or Global scale.

Conclusions
This study produced a soil erosion risk map for the most fire-prone vegetated areas in Portugal (eucalyptus, pine and shrubland), based on the revised Morgan-Morgan-Finney model predictions and validated using local post-fire field data, which, to the best of our knowledge, has never been done for Portugal.
The main conclusions drawn while developing this study can be summarized as follows: a. the highest post-fire soil erosion risk is estimated in north-central Portugal, and land managers should be aware of it, especially given the national wildfire recurrence history; b. the development of risk maps should consider different data sources to account for uncertainties in their predictions, thus creating better decision-making products; c. the validation of soil erosion estimations with field data is important to evaluate the robustness of model predictions, compare distinct scenarios, and understand which data sources are the most suitable for model application.