Modelling the Event-Based Hydrological Response of Mediterranean Forests to Prescribed Fire and Soil Mulching with Fern Using the Curve Number, Horton and USLE-Family (Universal Soil Loss Equation) Models

Abstract

The SCS-CN, Horton, and USLE-family models are widely used to predict and control runoff and erosion in forest ecosystems. However, in the literature there is no evidence of their use in Mediterranean forests subjected to prescribed fire and soil mulching. To fill this gap, this study evaluates the prediction capability for runoff and soil loss of the SCS-CN, Horton, MUSLE, and USLE-M models in three forests (pine, chestnut, and oak) in Southern Italy. The investigation was carried out at plot and event scales throughout one year, after a prescribed fire and post-fire soil mulching with fern. The SCS-CN and USLE-M models were accurate in predicting runoff volume and soil loss, respectively. In contrast, poor predictions of the modelled hydrological variables were provided by the models in unburned plots, and by the Horton and MUSLE models for all soil conditions. This inaccuracy may have been due to the fact that the runoff and erosion generation mechanisms were saturation-excess and rainsplash, while the Horton and MUSLE models better simulate infiltration-excess and overland flow processes, respectively. For the SCS-CN and USLE-M models, calibration was needed to obtain accurate predictions of surface runoff and soil loss; furthermore, different CNs and C factors must be input throughout the year to simulate the variability of the hydrological response of soil after fire. After calibration, two sets of CNs and C-factor values were suggested for applications of the SCS-CN and USLE-M models, after prescribed fire and fern mulching in Mediterranean forests. Once validated in a wider range of environmental contexts, these models may support land managers in controlling the hydrology of Mediterranean forests that are prone to wildfire risks.

1. Introduction

Wildfire is one of the most dangerous threats to forest ecosystems, since it impacts almost all components (air, soil, plants, fauna, surface water [1]). Fire alters the soil properties and removes vegetation, leaving the soil bare and, thus, exposed to flooding and erosion [2,3,4]. The hydrological impacts of fire are particularly important in Mediterranean forests. Here, the frequency, extent, and intensity of wildfires have been associated with an increase in climate warming in the last three decades [5], due to specific weather conditions (e.g., low humidity, high temperature, and strong winds [6]), and hydrological regimes (extreme and flash storm events with heavy and erosive rainfalls [7,8]).

To control fire severity and frequency, and at the same time mitigate the hydrological impacts in fire-affected areas, several pre- and post-fire management actions have been proposed. Prescribed fire applications, as well as soil mulching, have been used in several environments with satisfactory results [9,10]. Many studies have demonstrated how low-intensity fires are effective for removing the fuel that can generate a high-intensity fire, and therefore for reducing the wildfire risk in treated forests [11,12]. Soil mulching with vegetation residues is one of the most common post-fire management strategies to limit runoff and erosion in the short term, since mulch protects soil from rainsplash prior to vegetation regrowth [13,14].

The hydrological processes in burned soils are very complex, since several factors (weather, fire severity, vegetation cover, soil properties, morphology, and land management) influence the hydrological response of soil [2,4,15]. Computer-based hydrological models are essential tools to better understand and predict the hydrological processes in a cost-effective and time-efficient way [16], including in burned soils. However, validation under different environmental conditions is required, to ensure model applicability and reliability for predicting post-fire hydrology [17,18].

Many models have been developed to forecast runoff, erosion, and transport of polluting compounds under a large variety of climatic and geomorphological conditions [18,19]. Empirical models are of easier and quicker to apply compared to the most sophisticated physically-based routines. Their prediction capability may be considered acceptable for many needs and uses [20,21]. The use of empirical models is practically compulsory in data-poor environments, where input parameters are scarce or are difficult to collect, so that more complex models cannot be implemented [22,23]. The Soil Conservation Service (SCS)-curve number (CN) (henceforth “SCS-CN model”), Horton, and USLE-family (universal soil loss equation) models are the most common empirical methods adopted in various environments to predict runoff and erosion, respectively [18,19,24]. Thanks to the advantages of these empirical models, the SCS-CN, Horton, and USLE methods have been used as the hydrological and erosion components of many hydrological models applied at catchment scale (e.g., AnnAGNPS and SWAT models).

Regarding runoff and erosion predictions after fire using these models, several authors have evaluated the performance of hydrological models in burned forest areas, with or without post-fire treatments [4]. Regarding runoff prediction, the SCS-CN method is commonly used in fire-affected forests, but the CN values are still not well-known for burned conditions [25,26]. According to the latter authors, the literature on SCS-CN model accuracy and CN values for simulating runoff in burned forests is limited. Moreover, the variability of the hydrological predictions of different models, using the SCS-CN method as sub-component, highlights the difficulties in setting proper CN values for post-fire conditions using observed data [25].

The USLE-family models (USLE, RUSLE, MUSLE, USLE-M, HUSLE, etc.) have been developed to estimate soil erosion in agricultural lands, and their applicability in burned forest is not predictable, since fire impacts on vegetation and soil properties are different, and complicated by several factors (e.g., soil water repellency, ash effects), especially in the Mediterranean environment [27]. Moreover, these equations have mainly been applied for multi-year predictions, while modeling tests after the so-called ‘window of disturbance’ [28] (occurring immediately after a fire and throughout about one year, when the soil is bare and the pre-fire levels of soil properties have not yet been recovered) have been mainly carried out using RUSLE in soils burned by wildfires (e.g., [29,30]). To the authors’ best knowledge, no applications are available in the literature about the use of the MUSLE and USLE-M models after prescribed fires. As such, the optimal values of CNs and C factors, which are key parameters for the accurate estimation of runoff and erosion [31,32], have not yet been identified in burned conditions and with post-fire management. Moreover, the validation of soil erosion models in soils subjected to post-fire management treatments, such as mulching, is particularly scarce globally [28,33]. Therefore, the applicability of the SCS-CN method and USLE-family models may remain questionable, without targeted modeling evaluations.

These literature gaps require studies that should assess the prediction capability of the SCS-CN, Horton, and USLE models in burned forests in Mediterranean areas under pre-fire and post-fire management, such as prescribed fire and soil mulching.

To satisfy this need, this study evaluates the prediction capability for runoff and soil loss of the SCS-CN, Horton, and USLE-family models (MUSLE and USLE-M) in three forests (pine, chestnut, and oak) of Southern Italy. The investigation was carried out at the plot and event scales throughout one year after a prescribed fire and with post-fire soil mulching with fern residues. The research questions which this study aims to answer are two: (i) Are the tested models reliable and accurate for predicting surface runoff and soil erosion in Mediterranean burned forests? (ii) Which are the optimal values of the input parameters of the tested models?

2. Materials and Methods

2.1. Study Area

The investigation was carried out in three of the most dominant forests of Calabria (Southern Italy), whose climate is semi-arid (‘Csa’ class, ‘hot-summer Mediterranean’ climate, according to Koppen) [34]. The mean annual precipitation and temperature are 1102 mm and 17.4 °C, respectively (weather station of Sant’Agata del Bianco, geographical coordinates 42°17′54″ N, 59°51′59″ E, period 2000–2020).

Close to the municipality of Samo, three forest sites were identified to collect the hydrological observations used for the model evaluation (Figure 1):

• a pine (Pinus pinaster Aiton, “Calamacia” site, 38°4′52″ N; 16°1′46″ E) stand reforested in 1984 over an area between 650 and 700 m a.s.l.
• a natural oak (Quercus frainetto Ten., “Rungia” site, 38°5′20″ N; 16°0′39″ E) stand (900–950 m a.s.l.)
• a chestnut stand (Castanea sativa Mill., “Orgaro” site, 38°4′59″ N; 16°1′50″ E) about 30 years-old, between 700 and 750 m.

The tree density was about 950 (pine), 225 (oak), and 725 (chestnut) trees/ha. The tree height was 21 (pine), 10 (chestnut), and 18 (oak) m, while the breast diameter was 28, 20, and 41 cm, respectively. Shrub formations mainly consisted of Quercus ilex L., Rubus ulmifolius S., and Bellis perennis L. (pine forest); Cyclamen hederifolium and Bellis perennis L. (oak); and Rubus ulmifolius S., Pteridium aquilinum L., and Bellis perennis L. (chestnut). All forest stands had not been subject to management actions after planting or in the last fifty years for the natural stand.

The soils of the experimental sites (Cambisols, according to the World Reference Base for soil resources classification) were homogenous. The mean slope of soils was about 20% for all stands, and the texture was loamy sand (10.6 ± 2.57% of silt, 8.76 ± 0.61% of clay, and 80.7 ± 2.68% of sand). The unburned area of the pine forest in Calamacia instead showed a sandy loam texture (10.1 ± 1.01% of silt, 9.0 ± 0.01% of clay, and 81.0 ± 0.99% of sand).

2.2. Prescribed Fire Operations and Mulching Application

The prescribed fire was carried out in early June 2019, following the national and regional rules, and taking care that the main conditions during fire application at the experimental site (absence of wind and air humidity between 50 and 60%) were ideal, to avoid an uncontrolled wildfire. The burn severity of soils after the prescribed fire was low, according to the classification by [35].

Straw is commonly used as a mulch material, mainly in croplands, and this mulch is not always suitable in forest areas. Biomass transport from agricultural sites may be expensive, and these vegetal residues often contain agro-chemicals and parasites, with the possible development of non-native vegetation and diseases to forest plants [36]. Woodchips are sometime used as a mulch material in forests (particularly in young stands), but production may be expensive and difficult, due to the need of big machinery and a large amount of wood biomass [37]. Fern (Pteridium aquilinum (L.) Kuhn), which is abundant on the Mediterranean forest floor and does not bring non-native seeds or chemical residues into the forest ecosystem, may replace straw for the mulching of fire-affected areas. However, no studies prior to the present investigation have evaluated its suitability as a mulch cover for burned soils.

At the experimental sites, one day after the fire, small pieces (maximum length of 5 cm) of fern stems were applied to the soil as a mulch material in a part of the burned area. The fern was supplied from the same forest and the fresh residues were spread on the ground at a dose of 500 g/m² of fresh weight (200 g/m² of dry matter, as suggested for straw mulching by [38,39]).

2.3. Hydrological Monitoring

In each forest site, three series of plots (each one with three replications) were delimited at a reciprocal distance, between 1.5 and 20 m (Figure 2).

The plots (3 m in length × 1 m in width, for a total area of 3 m²) were hydraulically isolated, in order to prevent the inflow of surface water, using 0.3-m high metallic sheets inserted up to 0.2 m below the soil surface. Downstream of each plot, a transverse channel intercepted the water and sediment flows, which were collected in a 100-litre tank.

The hydrological monitoring started immediately after site installation and was carried out until September 2020, over 15 months (

Table 1. Main hydrological variables of rainfall events monitored at the experimental sites (Samo, Calabria, Southern Italy).

Date	Height (mm)	Net Height (mm) *			Duration (h)	Intensity (mm/h)
		Pine	Oak	Chestnut		Max	Mean
15 July 2019	65	61.8	59.8	60.5	36	22.2	1.99
9 October 2019	49.9	45.4	43.9	44.9	26	14.6	1.85
11 November 2019	142.8	135.7	132.8	132.8	41	26.2	3.49
23 November 2019	87.1	82.7	81.0	81.9	19	24.7	4.58
5 December 2019	147.2	141.3	138.4	139.8	30	19	4.90
24 March 2020	155.9	149.7	146.5	149.7	32	13.8	2.86
14 July 2020	22.4	20.6	19.7	20.4	7	12.8	2.58

Note: * recorded at the rain gauge station under a tree canopy in each forest.

). Precipitation height and intensity were measured in 15-min steps by a tipping bucket rain gauge at a maximum distance of 1 km from the experimental sites. Surface runoff and sediment concentration after the monitored events were measured according to Lucas-Borja et al. (2019b) [40]. In short, after mixing the water in the tank, three separate samples were collected for each rainfall–runoff event (total of 0.5 L). The samples were oven-dried at 105 °C for 24 h in the laboratory. Then, the dried sediments were weighed, and the weight was divided by the sample volume, to calculate the sediment concentration. The product of the latter by the runoff volume gave the soil loss.

2.4. Short Description of the Models

Some brief information about the tested models is provided below, while more details are available in the works by the cited authors.

2.4.1. SCS-CN Model

The Soil Conservation Service-curve number (SCS-CN) [41] was developed by the United States Department of Agriculture in the 1950s. This empirical model derives some assumptions from physically-based infiltration equations and requires only a few data points to estimate runoff for a given rainfall event.

The SCS-CN method assumes:

$$\frac{V}{P_n}=\frac{W}{S}$$

(1)

where V is the runoff volume, P_n is the net rainfall, W is the soil potential retention, and S is the maximum soil potential retention (all values are in mm).

V is calculated by the following equation:

$$V=\frac{P_n{}^2}{P_n+S}$$

(2)

where P_n is the difference between the rainfall depth P and the initial abstraction I_a (both in mm). The latter is the amount of rainfall retained in soil storage as interception, infiltration, and surface storage before runoff begins [42]. By convention, I_a is equal to the product of a coefficient λ (generally equal to 0.2) by S.

Therefore, V becomes:

$$V=\frac{{\left(P-0.2S\right)}^2}{P_n+S}$$

(3)

S is a function of the dimensionless ‘curve number’ (CN) parameter:

$$S=25.4\left(\frac{1000}{CN}-10\right)$$

(4)

CN describes the antecedent potential water retention of a soil [43]. Theoretically, CN varies between 0 and 100, but the usual values of CNs are in the range 40–98 [42].

The CN of agro-forest soils depends on the soil hydrological class, vegetal cover, hydrological condition (good, medium, poor), and cultivation practice; moreover, for CN calculation the antecedent moisture condition (AMC) of the soil must be determined. The soil hydrological class (A to D) is related to the soil’s capability to produce runoff, which in turn is due to the soil infiltration capability. The actual AMC of the soil subject to a rainfall/runoff event is estimated as a function of the total height of precipitation in the five days before the event in the two different conditions of crop dormancy or growing season. In this regard, three AMCs are identified:

• AMC_I: dry condition and minimum surface runoff
• AMC_II: average condition and surface runoff
• AMC_III: wet condition and maximum surface runoff.

The SCS-CN guidelines report tables to calculate the CN values for soils of a given hydrological class and condition, vegetal cover, cultivation practice, and average AMC (AMC_II). The values of CNs related to AMC_I (CN_I) or AMC_III (CN_III) can be calculated with the following equations:

$$C{N}_I=\frac{4.2C{N}_{II}}{10-0.058C{N}_{II}}$$

(5)

$$C{N}_{III}=\frac{23C{N}_{II}}{10+0.13C{N}_{II}}$$

(6)

2.4.2. Horton Equation

Horton’s method was formulated by Robert E. Horton in 1939 as an infiltration model to describe the physical process of infiltration in a quantitative manner.

The runoff rate q (in mm h⁻¹) at a given time t is given by:

q(t) = i(t) − f(t)

(7)

where i(t) and f(t) (both in mm h⁻¹) are the rainfall intensity and infiltration rate at time t, respectively.

The infiltration rate f(t) is calculated as:

f(t) = f_c + (f₀ − f_c) e^−kt

(8)

During a storm, f(t) generally declines from the maximum rate f₀ to the minimum value f_c through the parameter k. Equation (7) gives q(t) when i(t) exceeds f(t). The runoff volume is the integral of Equation (7), when q(t) is positive, between the start and the end of the runoff event.

2.4.3. MUSLE Equation

The ‘universal soil loss equation’ (USLE) was first established in the USA to model erosion in small agricultural catchments. USLE has a mathematical form that depends on six input parameters linked to climate, soil cover and properties, topography, and human activities; the six so-called “USLE-factors” (R, K, L, S, C, and P).

The USLE equation has been modified and updated over several versions and has been replaced by the revised USLE (RUSLE) [44,45]. Reference [46] developed a modified version, called MUSLE, which is the acronym modified USLE. The MUSLE model replaces the USLE rainfall factor (R) by a runoff factor, to consider the effect of flow on sediment transport. Therefore, the expression of the MUSLE equation has the following general form:

Y = a (Q’q_p)^b K L S C P

(9)

where Y is the soil loss (tons ha⁻¹) on a storm basis, Q is the runoff volume (m³), q_p is the peak flow rate (m³ s⁻¹), K is the soil erodibility factor (tons h MJ⁻¹ mm⁻¹), L and S are the slope-length and steepness factor, respectively, C is the cover management factor, and P is the conservation practice factor. The a and b coefficients are site-specific empirical factors for calculating the runoff factor.

2.4.4. USLE-M Equation

Kinnell and Risse (1998) [47] proposed the USLE-M model based on the hypothesis that the sediment concentration in the runoff is affected by the event rainfall erosivity index (R_e, [48] per unit quantity of rain (P_e, mm). According to the USLE-M, Y is calculated as:

Y = Q_R R_e K LS C P

(10)

where Q_R and R_e are the runoff coefficient and the erosivity index for the modelled event, respectively. The other factors of the USLE-M have the same meaning as the USLE and MUSLE equations, but the values of K and C factors are calculated using different expressions (see Section 2.4.3 and Section 2.4.4) [47].

2.5. Model Implementation in the Experimental Plots

2.5.1. SCS-CN Model

The sub-hourly precipitation records collected at the rain gauge stations were aggregated in daily values and supplied as input to the SCS-CN model. The AMC was derived according to the antecedent rainfall depths of each precipitation event. The soil hydrological group was identified using the data of the soil map of Calabria [49] and according to [50], who measured the hydraulic conductivity of the same sites.

The default values of CN were assumed, following the standard procedure by the USDA Soil Conservation Service [41] (

Table 2. Values of input parameters adopted to simulate surface runoff volumes and soil loss using the SCS, Horton, MUSLE, and USLE-M models applied in the experimental plots.

Model	Input Parameter	Measuring Unit		Soil Conditions
				Unburned			Burned				Burned and Mulched
				Default Model		Calibrated Model	Default Model		Calibrated Model		Default Model	Calibrated Model
Chestnut
SCS-CN	CN	-		46		43	70		80 *	45	50	65 *	32
	λ	-		0.2
Horton	m	mm h−1		33.65		-	30.51		-		37.61	-
	n	s−1		0.006		-	0.004		-		0.004	-
	r2	-		0.90		-	0.95		-		0.99	-
MUSLE	a	-		89.6
	b	-		0.56
	K-factor	tons h MJ−1 mm−1		0.03
	C-factor	-		0.009		1	0.390		1		0.011	1
	P-factor	-		1
USLE-M	Qr	-	max	0.17			0.34				0.10
			min	0.07			0.09				0.02
	Re-factor	MJ mm ha−1 h−1	max	69.5			84.8				54.6
			min	2.86			3.26				1.05
	KUM-factor	tons h MJ−1 mm−1		0.043			0.024				0.102
	CUM-factor	-		0.004		0.021	0.203		0.288 *	0.038	0.008	0.019 *	0.004
	P-factor	-		1
Oak
SCS-CN	CN	-		46	45		80		83 *	47	50	76 *	45
	λ	-		0.2
Horton	m	mm h−1		17.95	-		16.38		-		22.93	-
	n	s−1		0.007	-		0.005		-		0.005	-
	r2	-		0.67	-		0.95		-		0.90	-
MUSLE	a	-		89.6
	b	-		0.56
	K-factor	tons h MJ−1 mm−1		0.03
	C-factor	-		0.001	1		0.356		1		0.011	1
	P-factor	-		1
USLE-M	Qr	-	max	0.19			0.48				0.35
			min	0			0.14				0.07
	Re-factor	MJ mm ha−1 h−1	max	62.7			84.8				72
			min	0			24.9				8.16
	KUM-factor	tons h MJ−1 mm−1		0.052			0.023				0.034
	CUM-factor	-		0.001	0.020		0.356		0.104 *	0.056	0.011	0.056 *	0.045
	P-factor	-		1
Pine
SCS-CN	CN	-		39	40		90	79 *		41	55	65*	39
	λ	-		0.2
Horton	m	mm h−1		34.69	-		39.90	-			32.44	-
	n	s−1		0.003	-		0.00	-			0.004	-
	r2	-		0.95	-		0.98	-			0.94	-
MUSLE	a	-		89.6
	b	-		0.56
	K-factor	tons h MJ−1 mm−1		0.03
	C-factor	-		0.001	1		0.36	1			0.003	1
	P-factor	-		1
USLE-M	Qr	-	max	0.08			0.34				0.10
			min	0			0.06				0.06
	Re-factor	MJ mm ha−1 h−1	max	50.1			75.8				33.8
			min	0			15.6				10.9
	KUM-factor	tons h MJ−1 mm−1		0.085			0.033				0.068
	CUM-factor	-		0.004	0.008		0.203	0.208 *		0.006	0.008	0.055 *	0.003
	P-factor	-		1

Notes: CN = curve number; λ = initial abstraction ratio; m, n = coefficients of Equation (13); r² = coefficient of determination; a, b = site-specific factors; R_e = RUSLE rainfall R-factor; K, C and P = factors of the MUSLE model; K_UM and C_UM = factors of the USLE-M model; * first two modeled events.

).

2.5.2. Horton Equation

In the same experimental sites, [50] determined the water infiltration curves for the three soil conditions using a rainfall simulator (Eijkelkamp^®, https://en.eijkelkamp.com/), following the methods reported by [51]. In short, for each forest stand and soil condition, rainfall simulations were carried out in three randomly chosen points. Rainfall of 3.0 mm, at an intensity of 37.8 mm/h, was generated over a surface area of 0.305 m × 0.305 m. Throughout the simulated rainfall, the surface runoff volume was collected and measured in a small graduated bucket at a time scale of 30 s. The infiltration curves were determined by subtracting the runoff from the rainfall at each time interval. The infiltration test stopped when three equal time measurements of instantaneous infiltration had been recorded.

For Equation (8), we interpolated these infiltration curves using Equation (13), which has the following mathematical structure:

$$f\left(t\right)=m{e}^{-nt}$$

(11)

where m and n are the two constant coefficients and t is expressed in seconds. The goodness-of-fit of this equation was measured by the coefficient of determination (r²) (Table 2).

For the modeled events, the hyetograph i(t) was derived from the rainfall records and the difference between i(t) and f(t) at a given t gave the runoff rate q(t) every five minutes. Given the very short time of concentration (less than one minute) of the plot, the surface runoff stop was considered the same as the rainfall end.

2.5.3. MUSLE Equation

The MUSLE model is usually applied at the catchment scale; however, in some studies, it has been implemented at the plot scale (e.g., [52,53]). Q of Equation (9) was the runoff volume predicted by the SCS-CN method. The parameter q_p was calculated using the following formula:

$$q_p=z\frac{AQ}{t_p}$$

(12)

where A is the plot area (m²), z is a conversion coefficient, and t_p (0.01 h in this study) is the plot concentration time, which was experimentally measured at the plots using a surface tracer.

We deliberately adopted the values of Q and q_p as predicted by the SCS-CN method or calculated using Equation (12), rather than using the observed values, since these observations are never available in practical applications, and, therefore, the modeler is forced to use estimations from models. The site-specific factors, a and b, in Equation (8) were 0.87 and 0.56, according to the suggestions by [54].

The K-factor was estimated using the nomograph in [48]. The C-factor was calculated using an empirical equation based on canopy cover and aboveground biomass proposed by [55]. Usually, the C-factor is calculated as the product of some sub-factors, depending on the previous land use, canopy cover, surface cover, surface roughness, and soil moisture content (e.g., in [56]). In this study, we preferred to use a simpler equation, based on soil cover, since often these sub-factors are not easy to measure or determine. The P-factor was always set to one (Table 2).

2.5.4. USLE-M Equation

The runoff coefficient Q_R of Equation (10) was calculated as:

Q_R = Q/P_e

(13)

where Q and P_e are the runoff volume (again estimated using the SCS-CN method) and rainfall depth (both in mm), respectively. Following Renard et al. (1991) [45], the rainfall R-factor (R_e, MJ mm ha⁻¹ h⁻¹) factor was calculated using the following equation:

R_e = EI₃₀/1735

(14)

where NSE is the rainfall kinetic energy (tons m ha⁻¹) and I₃₀ is the maximum rainfall intensity in 30 min (mm h⁻¹) for each event.

The correction of K proposed by Kinnel et al. (1998) was adopted to calculate the K-factor of USLE-M (K_UM). The C- and P-factors (C_UM) were calculated as for the MUSLE equation (Table 2).

2.6. Model Calibration

The SCS-CN, MUSLE, and USLE-M models were initially run with default parameters, estimated from the guidelines of the three models. Since the models provided poor predictions (see Section 3), we chose the most sensitive input parameters (CN for the SCS-CN model [57,58], and the C-factor for the MUSLE and USLE-M equations [59,60]) and calibrated the models. In our study, the hydrological effects of prescribed fire and soil mulching were considered by tuning these parameters. The models were first calibrated using constant CNs and C-factors over time. In a further calibration trial, the CNs and C-factors were increased for the first two rainfall events, to take into account the variability of the soil hydrological response throughout the first year after fire, as demonstrated by several authors (e.g., [22,61,62]).

Calibration was carried out manually using a trial-and-error procedure, and was considered optimal when the coefficient of efficiency (NSE, see Section 2.6) was the highest and the error between the mean values of the observations and simulations of runoff (SCS-CN and Horton models) or soil loss (MUSLE and USLE-M models) was the lowest among the runs. The Horton equation was not calibrated, since all the relevant input parameters come from measurements from the infiltration tests. Table 2 reports the default and calibrated parameters, with their sources, for the four models.

2.7. Model Performance Evaluation

The runoff/erosion prediction capability of the four models was evaluated using qualitative and quantitative approaches. First, the observed and simulated values were visually compared around the line of perfect agreement in scatter plots. Then, we adopted a combination of the following criteria for model quantitative evaluation: (i) the main statistics (i.e., maximum, minimum, mean, and standard deviation of both observed and simulated values); and (ii) a set of indexes, commonly used in hydrological modeling. These indexes consisted of the determination coefficient (r²), the efficiency coefficient (NSE, [63]), and percentage bias (PBIAS, [64]), adopting a p-level of 0.05. In more detail, r² measures the dispersion of the ‘observations vs. predictions’ points around the interpolating line; values over 0.5 are deemed acceptable [56,65,66]. NSE, which is the ‘goodness of fit’ of the model predictions, is optimal if NSE = 1, good if NSE ≥ 0.75, satisfactory if 0.36 ≤ NSE ≤ 0.75, and unsatisfactory if NSE ≤ 0.36 [66]. PBIAS indicates whether the model over-predicts (if negative) or under-predicts (if positive) the output variable. A PBIAS below 0.25 and 0.55 for runoff and erosion, respectively, are considered fair [67,68].

3. Results and Discussions

3.1. Hydrological Characterization

Throughout the monitoring period, 516 rainfall events with a total depth of 1120 mm were recorded at the rain gauging stations. Of these events, only seven were classified as erosive events (that is, with a depth over 13 mm), according to [48] (Table 1).

In the unburned plots, the maximum runoff (up to 18.1 ± 12.9 mm in chestnut forest) was observed after the highest rainfall (156 mm). In two events (9 October 2019 and 14 July 2020), no runoff was collected in the unburned chestnut and oak forests and in all soil conditions of the pine and oak sites (Figure 3).

The highest runoff volume in the burned plots was always collected after the first post-fire event. One month after the prescribed fire, the runoff was from 22.3 ± 1.35 mm (chestnut) to 31.3 ± 2.29 (oak) (Figure 3).

The first rainfall event also produced the highest runoff in the burned and mulched plots of the chestnut and oak forests (6.61 ± 1.16 mm and 23 ± 3.69 mm, respectively). Conversely, in the pine forest, the maximum runoff (10.4 ± 0.80 mm) was measured after the second event (Figure 3).

The soil loss in the unburned plots was in the range 5.31 ± 1.40 g/m² (pine forest) to 15.34 ± 3.21 g/m² (chestnut), and these values were measured after the first and second post-fire rainfalls (Figure 4). For these two events, erosion increased very much in the burned soils of all forests, and particularly in the pine and chestnut soils (soil loss equal to 51.61 ± 6.92 and 52.26 ± 13.67 g/m², respectively). In oak soils, the soil loss was noticeably lower (15.12 ± 2.87 g/m²), but much higher compared to the unburned plots. Soil mulching with fern was effective for reducing erosion, and, under this soil condition, the maximum soil losses were between 10.62 ± 0.99 g/m² (pine forest) and 14.58 ± 4.80 g/m²; again recorded after the first event. After the first two events, the soil loss showed a low variability in the unburned soils, while, in the burned and not treated soils, erosion decreased over time (Figure 4).

Soil mulching with fern mainly reduced the erosion in pine and chestnut forests compared to the fire-affected plots. The maximum soil losses were equal to 1.87 ± 0.33 and 0.81 ± 0.16 g/m² (both surveyed in the third event), respectively. In these plots, the estimated soil losses were even lower compared to unburned soils, while the pre-fire erosion rates were only restored in oak forests for two events (Figure 4).

3.2. Hydrological Modeling

3.2.1. SCS-CN Model

The SCS-CN model, running with default input CNs, always gave poor predictions of surface runoff, as shown by the great scattering of the observations/simulations around the line of perfect agreement (Figure 5). This low accuracy is confirmed by the poor values of the evaluation indexes (

Table 3. Statistics and indexes evaluating the runoff prediction capability of the SCS-CN model in forest plots subject to prescribed fire and soil mulching with fern.

Run off Volume	Mean (mm)	Standard Deviation(mm)	Minimum (mm)	Maximum (mm)	r2	NSE	PBIAS
Pine
Unburned
Observed	7.00	4.54	0.00	13.06	-	-	-
Simulated (default)	5.18	5.74	0.00	12.34	0.73	0.36	0.26
Simulated (calibrated)	5.86	6.43	0.00	13.79	0.73	0.34	0.16
Burned
Observed	13.28	6.26	7.01	22.31	-	-	-
Simulated (default)	81.00	43.61	27.02	126.50	0.00	−85.36	−5.10
Simulated (calibrated)	12.41	7.06	0.52	22.28	0.69	0.81	0.07
Burned and mulched
Observed	7.41	2.31	4.37	10.35	-	-	-
Simulated (default)	20.04	18.32	0.32	40.58	−1.70	−81.26	0.79
Simulated (calibrated)	7.06	4.53	0.14	12.34	0.62	−0.70	0.05
Chestnut
Unburned
Observed	9.65	5.09	3.37	18.13	-	-	-
Simulated (default)	9.13	10.78	0.00	23.49	0.79	−0.72	0.05
Simulated (calibrated)	6.95	8.45	0.00	18.44	0.79	−0.13	0.28
Burned
Observed	15.39	6.39	3.85	22.31	-	-	-
Simulated (default)	35.35	31.85	0.00	74.04	0.14	−16.25	−1.30
Simulated (calibrated)	13.71	9.32	0.00	23.62	0.75	0.65	0.11
Burned and mulched
Observed	3.46	1.72	1.25	6.63	-	-	-
Simulated (default)	12.43	14.07	0.00	30.76	0.00	−5.14	−2.59
Simulated (calibrated)	2.86	2.77	0.00	8.13	0.72	0.94	0.17
Oak
Unburned
Observed	10.66	6.18	0.00	16.36	-	-	-
Simulated (default)	10.65	10.96	0.00	23.49	0.49	−0.66	0.00
Simulated (calibrated)	9.77	10.17	0.00	21.77	0.49	−0.43	0.08
Burned
Observed	20.18	7.09	10.00	31.34	-	-	-
Simulated (default)	59.25	37.76	13.74	99.21	0.04	−19.13	−1.94
Simulated (calibrated)	19.04	8.85	2.81	27.96	0.52	0.70	0.06
Burned and mulched
Observed	14.84	7.13	3.27	22.98	-	-	-
Simulated (default)	14.51	14.19	0.00	30.76	0.19	−1.28	0.02
Simulated (calibrated)	14.53	7.32	1.86	21.77	0.54	0.61	0.02

Notes: r² = coefficient of determination; NSE = coefficient of efficiency; PBIAS = coefficient of residual mass.

). In more detail, r² was much lower than 0.5 (with two exceptions, unburned soils in pine and chestnut forests, r² of 0.73 and 0.79), and NSE was below 0.35 (except for unburned soils in pine forest, NSE = 0.36). PBIAS, which was positive in some soil conditions and negative in others, indicates a high underprediction or overestimation for a observation, respectively. Moreover, the statistics calculated for the observations and predictions were highly different (mean error of up to 500%).

These model’s poor performances indicate that the literature values that were adopted as defaults for the input CNs were not suitable for simulating the runoff volume. This especially holds true under burned soil conditions, when the soil’s hydrological response was increased by fire. The CNs of fire-affected areas are usually estimated by increasing the post-fire values, depending on the fire severity (e.g., [69]). This statement agrees with the findings of [26], who highlighted the need to increase the CN values of burned soils by about 25 units. The worsening of the hydrological response of the burned soil after fires with different intensities has been shown by various studies (e.g., [70,71]), and this is also true in the case of prescribed fire (e.g., [1]). This increase is mainly due to the soil hydrophobicity and removal of vegetation due to fire. However, these effects vanish some months after a fire. The mulching treatment ‘smooths’ the increased hydrological response of burned soils, and this effect requires a lower increase in CN values [22].

However, in our experience, model runs with higher but constant CNs increase the runoff prediction capability of the SCS-CN model, but this calibration effort did not significantly improve the model performance in all soil conditions (data not shown). Regarding the unburned plots, the model predictions of runoff were also disappointing. Although the values of r² were satisfying (over 0.49), the NSE was negative in chestnut and oak forests, and lower than 0.36 in pine soils, while PBIAS (0.08) showed a noticeable tendency to model underprediction (Table 3). The difference between the mean values of the observed and predicted runoff was 8.3% to 28%. In contrast to our results and those by [26], references [27,72] observed no apparent increase of CNs in post-fire conditions.

In order to simulate the effects of a repellent and bare soil on surface runoff, there is a need to increase the CNs in the window-of-disturbance [3,28]. The authors of [22], working in burned pine forests under semi-arid Mediterranean conditions, demonstrated that the SCS-CN model performs better for simulating surface runoff, if the CNs are increased in the few months after the prescribed fire. Accordingly, when the CNs were increased for the first two modeled events after fire in our study, the performances of the calibrated model were always satisfactory for burned soils (mulched or not), except for the mulched soils in pine forest. The scattering of observations/simulations around the line of perfect agreement was reduced (Figure 5), and the evaluation indexes were generally over the acceptance limits for model predictions (Table 3). The SCS-CN model performance was good in burned soil for pine (r² = 0.69, NSE = 0.81 and PBIAS = 0.07), and satisfactory, both in burned and not treated, and burned and mulched soils, of chestnut (r² 0.72, NSE 0.65 and PBIAS 0.17) and oak (r² 0.52, NSE 0.61 and PBIAS 0.06) forests (Table 3). The difference between the mean observed and predicted runoff was lower than 17.5%. The model prediction capability of runoff was excellent in the burned and mulched soils of chestnut, where the r² and NSE were 0.72 and 0.94, respectively (Table 3).

The worse performance of the SCS-CN model in unburned soils compared to burned conditions was quite surprising. This low accuracy could be explained by the low generation capacity of the unburned soils of the experimental sites, which was not well simulated by the SCS-CN model, unless unrealistic CNs are input (runoff of three rainfall events simulated as zero) (Figure 5). The hydrological models generally tend to overestimate the lower events and underestimate the most intense flows [73,74]. A modified hydrological response due to fire increases the runoff generation capacity, which is better reproduced by this method. Overall, since the SCS-CN method does not optimally simulate the changes in soil properties due to management or other factors, without calibrating the input CNs, further studies should improve the model simulation of the temporal evolution of soil properties [75]. This could be done, for instance, by tuning the CNs proposed in the SCS guidelines using correction factors that should take into account the effects of soil water repellency and changes in hydraulic conductivity [76,77]. Until then, our results indicate that the suggested values of CN should be used instead of the standard SCS values for runoff predictions in soils burned by prescribed fires and treated with mulching under similar properties, climate, and management conditions as our experimental sites.

3.2.2. Horton Model

The runoff prediction capability of the Horton model was inaccurate under all soil conditions and forest species. In more detail, despite the satisfactory coefficients of determination calculated in the unburned soils of the three forest species (r² 0.65), the r² was always lower than 0.14 in the other soil conditions. The differences between the mean observed and predicted runoff volumes were over 50%, with peaks of up to 677%. Moreover, NSE and PBIAS were negative for all modeled soil conditions at the three sites (

Table 4. Statistics and indexes evaluating the runoff prediction capability of the Horton model in forest plots subject to prescribed fire and soil mulching with fern.

Run Off Volume	Mean (mm)	Standard Deviation(mm)	Minimum (mm)	Maximum (mm)	r2	NSE	PBIAS
Pine
Unburned
Observed	10.70	4.68	4.69	18.13	-	-	-
Simulated	24.83	17.05	6.12	41.63	0.65	−18.32	−1.32
Burned
Observed	17.31	4.24	12.49	22.31	-	-	-
Simulated	26.96	17.23	6.91	44.10	0.03	−5.43	−0.56
Burned and mulched
Observed	3.83	1.55	2.51	6.63	-	-	-
Simulated	26.51	18.49	1.71	44.48	0.14	−15.67	−5.92
Chestnut
Unburned
Observed	9.65	5.09	3.37	18.13	-	-	-
Simulated	25.03	18.30	1.92	45.49	0.76	−17.30	−1.59
Burned
Observed	15.39	6.39	3.85	22.31	-	-	-
Simulated	23.99	18.28	1.51	44.57	0.12	−3.81	−0.56
Burned and mulched
Observed	3.46	1.72	1.25	6.63	-	-	-
Simulated	24.06	17.92	0.87	44.21	0.04	−15.94	−5.95
Oak
Unburned
Observed	10.70	4.68	4.69	18.13	-	-	-
Simulated	31.51	15.43	14.95	46.50	0.68	−29.15	−1.95
Burned
Observed	17.31	4.24	12.49	22.31	-	-	-
Simulated	30.60	15.70	13.12	46.03	0.02	−6.06	−0.77
Burned and mulched
Observed	3.83	1.55	2.51	6.63	-	-	-
Simulated	29.77	16.17	11.30	45.80	0.08	−17.39	−6.77

Notes: r² = coefficient of determination; NSE = coefficient of efficiency; PBIAS = coefficient of residual mass.

). These coefficients indicate the wide scattering of the observations/simulations around the line of perfect agreement (Figure 6), and a noticeable overestimation of the modeled runoff volumes (Table 4).

The inaccuracy detected in this study for the Horton model was basically due to the noticeable overestimation of the observed runoff volumes, since the model was not calibrated, and, furthermore, the infiltration rate curves were not updated to account for the variability of infiltration rates over time. Furthermore, it cannot be excluded that the worse runoff prediction capability shown by the Horton model in comparison to the SCS-CN model may be due to the fact that the dominant runoff generation mechanism of the experimental soils is soil saturation during a storm (which is better simulated by the SCS-CN model) rather than infiltration excess (on which the Horton equation is based); as would be expected in soils in semi-arid environments [39].

3.2.3. MUSLE Model

Before calibration, the erosion predictions provided by the MUSLE equation were poor. The simulated soil losses were lower by at least one order of magnitude compared to the corresponding observations. The model strongly underestimated erosion in all soil conditions and forests (see the PBIAS close to one), which caused a great scattering of observations/simulations around the line of perfect agreement (Figure 7). These poor predictions were confirmed by the very low NSE and r² ( 0 and 0.20, the latter except for burned and mulched soils of chestnut forests, r² = 0.79) (

Table 5. Statistics and indexes evaluating the runoff prediction capability of MUSLE model in forest plots subject to prescribed fire and soil mulching with fern.

Soil Loss	Mean (tons/ha)	Standard Deviation(tons/ha)	Minimum (tons/ha)	Maximum (tons/ha)	r2	NSE	PBIAS
Pine
Unburned
Observed	0.02	0.02	0.00	0.05	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.04	−0.73	1.00
Simulated (calibrated)	0.01	0.01	0.00	0.01	0.04	−0.55	0.76
Burned
Observed	0.12	0.20	0.01	0.52	-	-	-
Simulated (default)	0.00	0.00	0.00	0.01	0.53	−0.06	0.97
Simulated (calibrated)	0.01	0.01	0.00	0.02	0.53	−0.01	0.92
Burned and mulched
Observed	0.02	0.04	0.00	0.11	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.01	−0.37	1.00
Simulated (calibrated)	0.01	0.00	0.00	0.01	0.01	−0.21	0.77
Chestnut
Unburned
Observed	0.05	0.05	0.00	0.15	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.18	−0.90	1.00
Simulated (calibrated)	0.00	0.00	0.00	0.01	0.18	−0.86	0.95
Burned
Observed	0.16	0.21	0.00	0.52	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.20	−0.25	0.99
Simulated (calibrated)	0.01	0.00	0.00	0.01	0.20	−0.22	0.97
Burned and mulched
Observed	0.03	0.05	0.00	0.15	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.79	−0.05	1.00
Simulated (calibrated)	0.00	0.00	0.00	0.00	0.79	−0.02	0.97
Oak
Unburned
Observed	0.05	0.02	0.00	0.07	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.05	−4.15	1.00
Simulated (calibrated)	0.01	0.01	0.00	0.02	0.05	−2.86	0.83
Burned
Observed	0.11	0.03	0.07	0.15	-	-	-
Simulated (default)	0.01	0.00	0.00	0.01	0.03	−1.35	0.94
Simulated (calibrated)	0.02	0.01	0.00	0.02	0.03	−1.01	0.86
Burned and mulched
Observed	0.07	0.03	0.04	0.12	-	-	-
Simulated (default)	0.00	0.00	0.00	0.00	0.00	−3.63	1.00
Simulated (calibrated)	0.01	0.01	0.00	0.02	0.00	−2.47	0.84

Notes: r² = coefficient of determination; NSE = coefficient of efficiency; PBIAS = coefficient of residual mass.

).

This inaccuracy of the MUSLE equation for simulating soil loss suggested the need to calibrate the model. We were forced to input C-factors equal to one, trying to reduce the model underestimation. This attempt was however disappointing, since the reliability of the calibrated MUSLE remained unsatisfactory. In addition, after calibration, the differences between predictions and observations were high, over 76% (Figure 7), and the evaluation indexes were poor. The values of NSE and r² were negative and lower than 0.20 in all forests (except r² = 0.79 in burned and mulched soils of chestnut), and the underestimation of soil loss was always high (as shown by the positive PBIAS, 0.76) (Table 5).

The unsatisfactory performance of the MUSLE model was quite surprising, despite calibration. We verified whether the applications of the Q and q_p simulated by the SCS-CN method rather than the observed values had influenced the erosion prediction capability of the MUSLE model. In modelling soil loss, the use of measured runoff volume and peak flow is suggested to improve the model accuracy, although this practice is often impossible in ungauged plots without runoff monitoring devices. However, the erosion prediction capability of the MUSLE model did not noticeably improve using the observed values of the hydrological variables in the runoff factor, since r² and NSE were lower than 0.50 and negative, respectively, and PBIAS was always over the acceptance limit of 0.55 for erosion prediction, as suggested in the literature.

The over-prediction of the MUSLE model is common in some studies carried out in different environments and soil conditions [78,79,80]. These authors reported that the low prediction capability could often be due to the fact that the model is applied in contexts that are different from the environments where the MUSLE was developed. More generally, the authors of [81,82] highlighted that small soil losses are usually over-predicted by USLE-family models.

3.2.4. USLE-M Model

As found for the MUSLE equation, the erosion predictions using the uncalibrated USLE-M were inaccurate, as visually shown in the relevant scatter plots (Figure 8). All the values of the evaluation indexes were unsatisfactory, since r² was lower than 0.41, NSE was negative, and PBIAS indicated strong model underprediction or overprediction (|PBIAS| 0.74; except for unburned, as well as burned and mulched, soils of pine, with PBIAS equal to 0.43–0.44, and therefore acceptable). Moreover, the differences between the mean or maximum values of predicted soil losses and the corresponding observations were always higher than 40% (with one exception, unburned soils of chestnut, 29%) (

Table 6. Statistics and indexes to evaluate the runoff prediction capability of the USLE-M model in forest plots subject to prescribed fire and soil mulching with fern.

Soil Loss	Mean (tons/ha)	Standard Deviation (tons/ha)	Minimum (tons/ha)	Maximum (tons/ha)	r2	NSE	PBIAS
Pine
Unburned
Observed	0.02	0.02	0.00	0.05	-	-	-
Simulated (default)	0.01	0.01	0.00	0.04	0.14	−0.93	0.43
Simulated (calibrated)	0.02	0.03	0.00	0.07	0.14	−2.33	−0.15
Burned
Observed	0.12	0.20	0.01	0.52	-	-	-
Simulated (default)	0.57	0.38	0.02	1.08	0.12	−5.84	−3.65
Simulated (calibrated)	0.12	0.19	0.00	0.42	0.87	0.90	−0.01
Burned and mulched
Observed	0.02	0.04	0.00	0.11	-	-	-
Simulated (default)	0.01	0.01	0.00	0.03	0.01	−0.20	0.44
Simulated (calibrated)	0.03	0.03	0.00	0.09	0.81	0.81	−0.08
Chestnut
Unburned
Observed	0.05	0.05	0.00	0.15	-	-	-
Simulated (default)	0.06	0.09	0.00	0.23	0.14	−3.97	−0.29
Simulated (calibrated)	0.05	0.07	0.00	0.18	0.14	−2.57	0.00
Burned
Observed	0.16	0.21	0.00	0.52	-	-	-
Simulated (default)	0.49	0.41	0.00	1.18	0.05	−4.06	−1.98
Simulated (calibrated)	0.16	0.21	0.00	0.53	0.95	0.96	0.01
Burned and mulched
Observed	0.03	0.05	0.00	0.15	-	-	-
Simulated (default)	0.05	0.05	0.00	0.12	0.41	0.25	−0.98
Simulated (calibrated)	0.03	0.05	0.00	0.13	0.86	0.88	−0.17
Oak
Unburned
Observed	0.05	0.02	0.00	0.07	-	-	-
Simulated (default)	0.01	0.01	0.00	0.02	0.05	−3.13	0.87
Simulated (calibrated)	0.05	0.06	0.00	0.15	0.05	−4.67	−0.02
Burned
Observed	0.11	0.03	0.07	0.15	-	-	-
Simulated (default)	0.45	0.29	0.07	0.93	0.07	−40.70	−3.26
Simulated (calibrated)	0.12	0.05	0.04	0.17	0.23	0.67	0.08
Burned and mulched
Observed	0.07	0.03	0.04	0.12	-	-	-
Simulated (default)	0.02	0.01	0.00	0.04	0.05	−2.00	0.74
Simulated (calibrated)	0.07	0.02	0.04	0.09	0.56	0.78	0.05

Notes: r² = coefficient of determination; NSE = coefficient of efficiency; PBIAS = coefficient of residual mass.

). Moreover, for this erosion model, we ascribed this poor performance to the tendency of hydrological models to overestimate and underestimate the lower and higher soil losses, respectively [18,73,81]. According to [18], the tendency for USLE-family models to overpredict low soil losses could be improved by incorporating an erosivity threshold in precipitation that must be exceeded before any sediment is generated.

The USLE-M model inaccuracy was removed thanks to calibration for the conditions of burned, and burned and mulched soils of all the forest species, while the erosion predictions provided by the calibrated USLE-M equation were still unsatisfactory for the unburned plots. For the latter soil condition, r² was lower than 0.14 and the NSE was negative (Table 6). In contrast, these evaluation indexes were over 0.56 (r²) (except in the burned soil of oak, r² = 0.23) and 0.67 (NSE) in burned soils (mulched or not) of all forests, and the |PBIAS| was lower than 0.17. The latter index reveals that in some soil conditions and forest species the model generally underpredicted erosion (burned soils, treated or not, of oak, and burned plots of chestnut), while, in the other cases, a slight tendency for the overestimation of soil loss was found). Moreover, the values of PBIAS were well below the acceptance limit of 0.55 stated in the literature ([67,68], see also Section 2.6). In addition, for burned soils of oak, the erosion prediction capability of the USLE-M equation can be considered as satisfactory, although the r² was low (0.23). As a matter of fact, both the NSE and PBIAS indexes complied with the acceptance limits (NSE 0.36 and PBIAS 0.55), and the differences between the mean or maximum values of the observations and predictions was only 8.5%. This statement is a proof that sometimes r² may be misleading in model evaluation [64,83], since it measures the scattering of values around the regression line and not around the line of perfect agreement.

The contrasting performances of the USLE-M model in predicting erosion between unburned and burned soils contrasts with the findings of [84], who reported insignificant impacts on erosion estimates between burned and non-burned forests.

Overall, for the USLE-family models, a calibration process has been considered necessary by several authors for improving their prediction accuracy. For instance, [85] and [86], applying the USLE-M in plots in Western Sicily (Italy) and under different soil conditions, highlighted the importance of the calibration process to allow its adaption to the different climatic and edaphic conditions.

Since five (K, L, S, C, and P) of the six USLE-factors are common in the two models under each soil condition, it is possible to compare the effects of the R-factor on the predicted soil losses. This indicates that, under the experimental conditions, the soil loss is mainly due to rainsplash erosion. This fact derives from the better performance of the USLE-M model, whose rainfall erosivity is based on the R-factor, compared to the MUSLE model. In contrast, the lower accuracy in simulating erosion shown by the latter equation, which includes parameters related to surface runoff, indicates the minor role of the soil loss produced in our plots by overland flow, which determines particle detachment. However, these statements depend strictly on the small scale of the experimental plots (only a few square meters), as well as the low number of plots, and thus must be verified at larger scales. As carried out for the MUSLE model, we ran the USLE-M model using the observed Q_R rather than the value calculated using the runoff volume predicted by the SCS-CN model, but the erosion prediction capability of USLE-M did not appreciably improve (data not shown). This means that the model can be applied to ungauged plots (that is, without equipment to measure runoff and peak flow) without losing accuracy when predicting soil losses in burned soils.

As previously mentioned, the MUSLE and USLE-M models had not been applied to predict erosion in fire-affected areas (neither by wildfires nor low-severity fires) before this study. Therefore, the results of the present study cannot be compared to similar experiments in the same or other environmental contexts. Comparisons with published studies that have evaluated the erosion prediction capability of USLE-family models are thus possible with reference to the RUSLE model. The authors of [56] found that the latter equation is ideal for fast and simple applications (i.e., prioritization of areas-at-risk) in zones burned by wildfires in Portuguese pine forests; since, in this application, an NSE between 0.63 and 0.70, and PBIAS in the range −12% to −20% were achieved. In contrast, in burned forests of pine and eucalyptus, and in lands with post-fire management in Galicia (NW Spain), references [29,30] and [87] found poor predictions of soil loss. Negative or very low positive NSE and high PBIAS were detected by these authors, since the RUSLE model overestimated the soil losses. These authors attributed the poor performance of the model to three factors: (i) the use of an inadequate kinetic energy equation of rainfall for the Mediterranean climate; (ii) the fact that soil water repellency, a key factor in post-fire soil hydrology [15,88], is not explicitly considered in the RUSLE model; and (iii) the poor ability of C- and K-factors to reflect changes in soil properties induced by fire. The attempts at calibration by tuning the R- and C-factors [30] and by introducing a reduction factor for soil erodibility (to take into account stone cover of the experimental soils) [29], or accounting for high contents of soil organic matter [87], did not noticeably improved the RUSLE accuracy for predicting soil loss. These statements agree with those reported by [18], who achieved a low reliability with the RUSLE equation in modeling sediment yields throughout the first year after fires of different severity in Colorado (USA). The same authors reported that the model accuracy was not improved by increasing the K-factor of RUSLE (which should simulate a decrease in aggregate stability after fire). Moreover, they advised that the use of RUSLE in unburned forests may be troublesome (as was found in our study, given the basically poor performance of the MUSLE and USLE-M models in the unburned plots), because overland flow is not common [18,89]. Since the main effects of fire are the changes in soil properties and surface covers [2,12,90], which may induce noticeable changes in the K- and C-factors, all these studies concluded that these USLE-factors do not adequately describe soil modifications after fire. Our study has instead demonstrated that satisfactory predictions of erosion in soils burned by prescribed fire (with or without mulch treatment) can be provided by the USLE-M model when using a simple calibration of the C-factor and regardless of the specific process (e.g., soil water repellency, decrease in aggregate stability, depletion of organic matter content) that affects the post-fire soil hydrology. A limitation of this procedure is the need to calibrate the C-factor in each environmental context (that is, in soils of different types and textures) before applying the USLE-M model. However, this requirement seems to be an easy task, considering the low time and money requirements of a small plot installation, as well as the low monitoring period (less than one year), as in the present study. This consideration is in close accordance with [56] and [84], who suggested that estimations of C-factors for burned areas should be determined for each context and fire type, although previous works reported several indicative values [30,91]. Therefore, the estimation of the site-specific C-factors using locally measured data, as was the aim of our study, increases the erosion prediction accuracy of USLE-family models [92].

A limitation of our study is that the C-factors proposed in this investigation for the modeled soil conditions were calibrated using observations collected at only one study area. Furthermore, the low number of plots (due to budget limitations) in this study may explain the low accuracy of the tested methods for modeling runoff and soil losses, especially in the unburned plots. In some cases, this procedure can be misleading and must be validated in other environmental conditions or supported by external parameters. Bearing in mind the limitations of our study, the C factors proposed for the experimental soil conditions could be reliable, at least for burned soils of Mediterranean forests (treated with mulching or not), and fill the lack of similar values for USLE-M applications in the literature.

Overall, the investigation enriched the literature about post-fire hydrological modeling, which is not homogeneously distributed worldwide and is still far from being exhaustive [93,94]. The results of our study confirm the applicability of two very common hydrological models (SCS-CN and USLE-M). Moreover, the unsuitability of two other methods (Horton equation and MUSLE) was demonstrated, at least in various forests in Mediterranean areas, whose intrinsic conditions (e.g., very shallow soils, strong soil water repellency, and peculiar hydrologic regime) often make the available hydrological models unsuitable. These models were developed in other climatic contexts and not in fire-affected areas, and therefore may find limited applicability without targeted modifications [56]. This study supports the efforts of modelers in choosing the most suitable hydrological models after prescribed fires and mulching treatments in Mediterranean forests that are prone to risk wildfire. Moreover, our modeling experiment proposed important input parameters (CNs and C-factors) for loamy sand soils supporting broadleaf (oak and chestnut) and conifer (pine) forest species, and the proposed values can be used under the evaluated soil and climatic conditions.

4. Conclusions

This study has demonstrated the feasibility of the SCS-CN and USLE-M models for predicting surface runoff and erosion, respectively, in plots burned by prescribed fire and mulched with fern residues in three Mediterranean forests of pine, chestnut, and oak. In contrast, poor predictions of the modelled hydrological variables were provided by the models in unburned plots, and by the Horton and MUSLE models for all the soil conditions. This result answers the first research question about the feasibility of the four models for hydrological applications in Mediterranean burned forests.

Regarding the second research question about the optimal values of the input parameters for the tested models, the study has proved that a calibration process is a prerequisite for the tested models for accurate predictions of surface runoff and soil loss. After calibration, we suggest the following two sets of CNs and C-factors for applications of the SCS-CN and USLE-M models after prescribed fire and fern mulching in Mediterranean forests:

• for runoff modelling CNs of 80, 83, and 79 (burned soils) and 65, 76, and 65 (burned and mulched soils) throughout the two to three months after fire, and 45, 47, and 41 (burned soils) and 32, 45, and 39 (burned and mulched soils) in the following period, for chestnut, oak, and pine forests, respectively;
• for soil loss predictions C-factors of 0.521, 0.085, and 0.339 (burned soils) and 0.069, 0.060, and 0.040 (burned and mulched soils) immediately after a fire, and 0.052, 0.059, and 0.008 (burned soils) and 0.007, 0.050, and 0.001 (burned and mulched soils) in the following period, for chestnut, oak, and pine forests, respectively.

The different CNs and C-factors consider the variability of the hydrological response of soil after fire.

Despite the satisfactory accuracy shown by the SCS-CN and USLE-M models in this modeling experiment, these results are limited to the experimental conditions; however, they are encouraging for further applications of these conceptually simple and widely used models in analogous climatic and geomorphologic conditions. Further modeling studies should also enlarge the spatial scale from plots to watersheds. Once validated in a wider range of environmental contexts, these models may support land managers in controlling runoff and erosion in forest soils that are prone to hydrogeological risks.