Rill Erosion in Unpaved and Rock-Paved Roads after Wildfire in a Mediterranean Forest

Forest roads are often subject to intense runoff and erosion, and the rates can be increased by other disturbance factors, such as wildfires. Since scarce literature exists on the effects of wildfires on rill erosion of forest roads, this study presents the first results of a wider research, evaluating rill erosion in four different types of roads on a forest in Hellìn (Castilla-La Mancha, Central-Eastern Spain): unpaved roads made of native materials (soil found at the study site) and rock-paved roads, both built in unburned areas as well as unpaved and rock-paved roads, in fire-affected areas. In general, the unpaved roads are more subject to rill erosion compared to the rock-paved roads. In particular, the road of burned areas shows an erodibility that is higher by more than 200% compared to the unpaved and unburned roads, and even by about 400% compared to rock-paved roads (in both burned and unburned areas). A modeling approach based on distance linear models and distance-based redundancy analysis has identified the slope of road surface and upstream hillslope as well as the percent bare soil over the road surface as important input variables to predict rill erosion in future modeling experiences. All these variables may be easily measured by quick field surveys. Although the analytical approach of this study is limited to the geometric characteristics of erosion features, the results and the methods developed are useful to support the activity of land managers to better understand the magnitude of road erosion and to develop efficient measures for its control


Introduction
Roads are fundamental infrastructures for several activities related to forest management (e.g., timber harvesting, site preparation, fire management, insect and disease control), and recreational opportunities [1,2]. However, roads are an element of hydrological disturbances in forested watersheds, because these infrastructures alter sediment connectivity and generate water and sediment flows [3,4]. Forest roads generally are the dominant sediment sources from forestlands [5,6] and this makes soil erosion particularly important in forests, because natural erosion rates are usually very low [7,8].
Road construction creates several cutslopes and fillslopes, as well as ditches and culverts, which become features of erosion and deposition [9]. The alteration of natural profiles of forests hillslopes and surface and subsurface water flows, the low plant cover, and the compaction of soil on the roadbed [10] are the factors that explain these processes [11]. The erosion features of forest roads are sources of increased erosion with enhanced land degradation [5], since soil is compacted by traffic and vegetal cover is low or absent [12,13]. Soil compaction and scarce vegetation reduce water infiltration, which on its turn enhances infiltration-excess overland flow, even for low-intensity rain events [14,15]. The impacts of roads are directly related to runoff and erosion rates, due to the highly compacted road surface, even though roads usually represent a small proportion of forest areas [16,17].
Road erosion depends on surface characteristics, such as soil texture, ground cover, time since construction and maintenance activities [7,18] and this makes the evaluation of the erosion rates a complex task. Erosion on roads is mainly in rills [2], and the presence of rills is of prime concern for forest management operation, because this affects the driveability of a road and water and sediment flows downstream can adversely affect water quality and aquatic habitats [2,19]. As demonstrated by [20], the mean annual sediment production in roads is strongly related to rill density, since a high rill density may determine nearly 20 times more sediment than roads without any rills.
The erosion rates in forests with roads can be enhanced by other disturbance factors, such as wildfires [21]. Wildland fires cause extreme changes in the landscape that can drastically influence soil cover and properties, surface runoff and soil loss [22]. The effects of post-fire water and sediment flows often require that many structures must be treated following fires, because roads in areas burned at high severity will produce much more runoff and erosion compared to other roads [2]. As reported by [16], wildfires and road construction may drastically reduce infiltration rates, thus increasing surface erosion rates by several orders of magnitude. Therefore, evaluating the erodibility of forest roads has an increasing importance [12], especially in forests burned by wildfires. The existing literature has mainly estimated sediment delivery rates as a function of upslope erosion rates but failed to isolate forest roads as site of erosion [5]. In other words, while it is well known that unpaved roads determine less erosion compared to rock-covered roads, the soil erodibility has rarely been evaluated in roads built in burned forests. To the best of our knowledge, no studies have directly examined how wildfires affect rill erosion on forest roads [20]. Only [2] have assessed the interactions between fire severity and road segment characteristics affect rill erosion and sediment deposition in forests of Colorado (USA), but their analysis did not evaluate how much rill erosion rates differ by road age and surface characteristics. Similar studies are completely absent in Western Europe, where the semi-arid climate may increase the wildfire risk and forests are threatened by several degradation factors besides fire (e.g., [23]), such as erosion rates beyond the tolerance limit (e.g., [24]). Moreover, it is also interesting to evaluate in burned forests whether the road age and surface characteristics influence the erosion rates, but no relevant studies are available at our knowledge.
To fill these gaps, this study presents the first results of a wider research targeted to quantify rill erosion and deposition rates in roads of different types subject to wildfire compared to unburned roads. In more detail, rill erosion is evaluated in two different types of roads (unpaved roads built by native materials and rock-paved roads) after a wildfire and in other roads with the same characteristics but selected in unburned areas of a forest of Hellín (Castilla-La Mancha, Central-Eastern Spain). The analytical approach is limited to the geometric characteristics of erosion features, such as the rills, which are simpler to be evaluated compared to other factors on which erosion depends (e.g., soil properties, lithology, etc.). We hypothesize that roads made of unpaved materials are more prone to rill erosion and that, for roads with the same surface characteristics, the infrastructures in unburned areas are less erodible compared to the fire-affected infrastructures. This fact may be modulated by wildfire effect on paved and unpaved roads. Moreover, a first attempt to propose models for predicting rill erosion from geometric characteristics of roads is proposed.

Study Area
The mountain of Sierra de Los Donceles (Hellín, Province of Albacete, Castilla-La Mancha Region, Spain) is located in the southeast of the Iberian Peninsula (Figure 1a,b). This area is lodged between the Mundo and Segura Rivers to the north and south, respectively. Altitude ranges from 304 m to 808 m at the Sierra de Los Donceles line. The mountain lies among pre-Baetic mountain chains with limestone and dolomite outcrops alternating with marly intercalations that date back to the quaternary. According to the Soil Taxonomy System, soils belong to the order of Aridisols, and suborder of Calcids, and show a loamy and sandy loam soil texture. The climate is Mediterranean semiarid, and the area is located on the meso-Mediterranean bioclimatic belt [25]. The mean annual temperature and precipitation are 16.6 • C and 321 mm, respectively. The maximum precipitation is usually in October (44.5 mm) and another peak is in May (39.6 mm) (1990-2014, data by the Spanish Meteorological Agency AEMET). The dry period goes from June to September, with an air relative humidity below 50%. Vegetation consists of Quercus cocciferae-Pinus halepensis S. series, where Aleppo pine comprises most of the tree cover strata and kermes oak mostly occupies the shrub strata. From records of forest wildfires in Spain starting in 1968, the most recent wildfire in Sierra de Los Donceles occurred in July 2012, when the area was devastated with roughly 6500 ha of burned forest. In May 2020, we selected four roads ( Figure 2) in the area completely affected by the wildfire of 2012 to carry out the present study. Burn severity was classified as high severity, according to [26], for the entire study area.

Road Construction and Characteristics
Four roads were identified to carry out this study ( Figure 2): • (a) and (b) an unpaved and a paved road in an unburned area (hereinafter indicated with "U/u" and "U/r", respectively, where "U" stands for "unburned", "u" stands for "unpaved" and "r" stands for "rock-paved"); • (c) and (d) an unpaved and a rock-paved road in an area affected by the wildfire ("B/u" and "B/r", respectively, where "B" stands for "burned").
In more detail, the U/u and B/u roads were made of native materials (soils found at the study area), while the U/r and B/r roads were paved with rock materials (mainly gravel). These roads showed the same characteristics in terms of soil, geology, geomorphology and level of traffic, and the same effects of disturbance due to the wildlife of July 2012.
According to the Forest Service of Castilla-La Mancha Region, the unburned roads (U/r and U/r) were built in 1960, open to traffic for about 45 years, and maintained every 7 years, approximately. Work maintenance consisted of simple re-grading, with displaced gravel being pushed back into shape. The burned roads (B/u and B/r) were built in September 2013, ripping the road surface to a depth of approximately 0.4 m with a tracked bulldozer pulling three unwinged ripping teeth. In the case of rock-paved roads, construction started with the subgrade layer. Then, a 5-cm layer of gravel was homogenously distributed on this surface creating a crown in the center point of the road; a 4-5% gradation from the center to the edge of the road was not exceeded. The gravel used was composed by crushed stone, sand, and fines. Fines were silt or clay particles smaller than 0.075 mm, which acted as a binder. In burned areas, paved and unpaved roads have been not maintained until to date.

Experimental Design
The roads and rills over their surface were mapped in May 2020, about 8 years after fire. The surveyed length varied from 517 m in B/r roads to 547 m in U/r roads and from 455 m in B/u roads to 474 m in U/u road. According to [2,18], roads were divided into hydrologically distinct segments, each one being typically defined by the road prism including the road surface, cutslope and fillslope (if present). The hillslope draining onto the road should be considered into the segment, but this source of runoff and erosion is commonly ignored because of the high infiltration rates of forested areas [2].

Measurement of Geometric Properties and Soil Covers of Road Segments
The following geometric characteristics were measured for each type of roads using a common measuring tape (for the horizontal distance) and a level (for the vertical distance):  Table 1).  Then, the areal covers of vegetation (VC), rock (RC) and bare soil (BS) were measured using touch lengths in one transect each three meters of segment length. These covers were expressed in percentage over the total area of the road (Table 1).
In the segments of the four road types all rills were mapped. The rill geometry, the width and depth of each rill were measured, and the cross-sectional area was calculated by assuming a triangular shape. A rill was identified as such, when its depth and width were at least 10 mm. On each road the following geometric characteristics of rills were measured using the same methods as for roads: (1) Rill Length (RL, m); (2) Rill Average Depth (RAD, cm); (3) Rill Average Width (RAW, cm); (4) Rill Average Slope (RAS, %). These geometric data allowed estimating the volume of rill erosion (RE). Rill drainage density (RDL), the ratio between the total length of rills and the road length, was calculated for each road type. The delivery of road sediment depends on the hydrologic connectivity (C), which is defined as the linkage or connection between a runoff source and the receiving water. According to [2], the following connectivity classes were adopted: class 1-No drainage feature (negligible potential for sediment delivery); class 2-Drainage features less than 10 m long that did not extend to within 5 m of stream (very low potential for sediment delivery); class 3-Drainage feature more than 10 m long but does not extend to within 5 m of an ephemeral or permanent stream channel (low to moderate potential for sediment delivery); class 4-Drainage feature extends to within 5 m of the stream (the associated road segment is connected and very likely to deliver at least some runoff and sediment to the channel network).

Statistical Analyses
First, the statistical differences in rill erosion volumes were determined by the multivariate permutational analysis of variance (PERMANOVA, [27]), using the four road types as factors. PERMANOVA tests the simultaneous response of one variable to one or more factors in an experimental design based on any resemblance measure, using permutation method. Before PERMANOVA, the geometric characteristics of roads and soil covers (the so-called "environmental variables") were log(x + 1) transformed, whereas RE volumes were square root transformed; these transformations are suggested by the PRIMER manual (the software package used for the statistical analyses), since transformations are usually applied as a method of changing the relative emphasis of the analysis on scarce versus more abundant data among categories.
The resemblance matrix was built using the Euclidean and Bray Curtis distance for environmental and RE data, respectively. The sums of squares type were type III (partial) and the four-level factor was a fixed effect (road type). The permutation method used was the unrestricted permutation of raw data and the number of permutations was 999.
Secondly, the non-metric Multi-Dimensional Scaling (NMDS) and the Kruskal stress formula (minimum stress: 0.01) were applied to the environmental variables and RE data, to evaluate the similarity level in the individual cases of the dataset.
Thirdly, a DISTLM function (DISTance-based Linear Modelling) was developed, to determine the relative importance of road characteristics and soil covers on RE volumes. For the DISTLM routine, we developed "marginal" tests of the relationship between the response variable (RE volumes) and an environmental variable (road characteristics and soil cover) to identify the independent variables that explain the variations in the erosion/deposition volumes. Following the marginal tests, "sequential" tests of individual variables were performed to assess whether adding an individual variable contributes significantly to the explained variation of the response variable.
Fourthly and finally, the distance-based redundancy analysis (dbRDA) was applied to RE volumes, to build a regression model against two new response variables ("axis" 1 and "axis" 2), built on the environmental variables. The AICc (Akaike Information Criterion, [28]) criterion was adopted to select the best model and the stepwise procedure was followed to build the model. These techniques were preferred to the traditional statistical tools used in other similar investigations (e.g., ANOVA, PCA and multi-regression models), since the latter, although being generally robust and meaningful, require some constraints in the input data (e.g., normality of data distribution or independent samples). The non-parametric techniques are lower limitations but are less powerful and significant. Conversely, PERMANOVA, NMDS, DISTLM and bRDA, based on the distance of data samples, can overcome these limitations of the traditional statistical techniques, but at the same time are robust and meaningful [29]. For the statistical analyses, the software PRIMER V7 ® with PERMANOVA add-on [27] and Statgraphics Centurion XVI ® (StatPoint Technologies, Inc., Warrenton, VA, USA) were used. A significance level of 0.05 was used, unless otherwise indicated.

Results and Discussion
The volume of eroded soil per unit of road length was the highest for B/u roads (0.0043 ± 0.0076 m 3 /m) and the lowest for U/r segments (0.0007 ± 0.001 m 3 /m) (Figure 3). According to PERMANOVA, rill erosion was significantly different (Pseudo-F = 3.39, p < 0.05) among groups. Moreover, the pairwise comparisons showed significant differences between B/u and U/u roads (t = 2.69, p(perm) < 0.001), U/r and U/u roads (t = 2.16, p(perm) < 0.021) and U/u and B/r roads (t = 1.87, p(perm) < 0.031). Instead, the comparisons among the other road types were always not significant (p(perm) > 0.05). This means that, after wildfire, the burned roads made of unpaved materials are much more erodible compared to the other road types, and rill erosion may increase by more than 200% compared to rock-paved, and unburned and unpaved roads, and even by about 400% compared to only rock-paved roads. As demonstrated by different authors, the unpaved forest roads significantly alter hillslope hydrology by increasing and concentrating surface runoff and erosion [16,30]. Moreover, the unburned roads, whose construction is older compared to the more recent infrastructures built in burned areas, show less erodibility, since the compaction level of the older roads is higher due to traffic and weathering, and this may have led to a lower sediment detachability from their surface ( Figure 3).
Our data show that after high and moderate severity fires the rills and gullies generated over the road surface extend for many tens or even hundreds of meters. The much longer lengths of these drainage features are due to the increased amounts of runoff and sediment flows from upslope, the accumulation and discharge of this runoff at a single location due to the road segment, and the reduced infiltration and roughness of the burned hillslopes below the road. Unpaved roads are recognized to have a heavy effect on several hydro-geomorphic processes, acting as a primary sediment source [7,12]. This is due to the lack of surface cover (differently from rock-paved roads), which exposes the road surface to both rainsplash, and overland and concentrated erosion [21].
The highest density of rill drainage network was found in U/r roads (0.53 ± 0.51 m/m), about 80% more compared to the lowest RDL, detected in B/r roads (0.27 ± 0.33 m/m), while B/u and U/u roads had intermediate values (0.47 ± 0.37 and 0.48 ± 0.44 m/m, respectively) (data not shown). The highest erosion detected in U/r roads (having the lowest RDL) is due to the largest cross-section, and, therefore, although the network is linearly developed, the sediment mobilised is low. In general, the more erodible unpaved roads showed a mean RDL that was 20% higher compared to the rock-paved roads. This larger erodibility of the surface in unpaved roads should be paid attention compared to the other road components (cutslopes and fillslopes), since the road surface generally produce the dominant source of runoff and sediment at each site [31].
NDMS provided two derivative components (here called MDS1 and MDS2), which explain together 63.2% of the total variance of the original variables (road geometry and soil cover) (Figure 4). In more detail, all soil covers have high loadings on MDS1 (0.879 for R/BS, −0.851 for R/RC and 0.807 for R/VC). R/US (loading of −0.545), C (−0.731), R/US (0.535) and R/DS (−0.369) have instead significant loadings on MDS2. This statistical analysis therefore helps to identify two new meaningful variables, of which one (the MDS1) is clearly related to soil cover and the other (MDS2) to road geometry, and both discriminate the different characteristics of the four road types. As a matter of fact, the same NMDS shows that two evident clusters of road segments can be identified, which group unpaved (burned or unburned) roads. Conversely, the segments of the other two road types are more scattered, with a partial overlay among several segments of U/r on both B/r and U/u roads ( Figure 5), and therefore no clusters can be evidenced. Another important consideration is that in some road segments the geometric characteristics and soil covers of the rock-paved roads may be similar as those of unpaved roads and this may indicate a higher erodibility.
By applying the distance linear models (DISTLM), the marginal tests revealed that only ROR/S has a significant (p < 0.001) influence on RE, while the other variables are much less significant (for R/DS, p < 0.1) or not significant (p > 0.1) at all, when those variables were considered as isolated ( Table 2). The sequential tests of DISTLM showed that the best distance linear model (R 2 = 0.76; AICs = 96.5) for predicting RE consisted of R/US and R/DS beside ROR/S ( Table 2). These results indicate that the road slope is the main driver of the rill erosion process, since the water flow energy eroding the road surface is primarily a function of the slope of the free surface stream. The main role of road slope in influencing rill erosion detected in this study is in close accordance with findings of other authors; for instance, [7,18,32] showed that road slope (often in combination with other geometric variables) can be used to predict road surface erosion. Moreover, also the upstream slope and the bare soil percentage have a significant influence on the concentrated erosion. Higher slopes increase the erosion process in upstream areas, while the road surface, if bare, is exposed to rainsplash and overland erosion. Additionally, these results are in line with other authors: [6,9] showed that slope properties of the roadsides as well as the conditions of road surface play a heavy influence on runoff and sediment flow generation in roads. Conversely, the other factors, such as the vegetation and rock cover as well as connectivity, play a much lower and not significant role in driving rill erosion of the studied roads. This may be quite surprising, since some researchers have demonstrated that both plant cover and presence of rocks have significant effects on runoff and sediment yields of roads (e.g., [33,34]). This apparent contrast may be explained by the fact that the roads analysed in this study have been subjected to a wildfire, which has completely removed the vegetal cover on their surface, and the differences in rock fragment covers among the four road types are too small for playing an influence on rill erodibility.
According to the variations (out of the fitted model and out of the total variation) explained by the axes of dbRDA, axis one (dbRDA1) applied to RE explained 99.32% of the fitted model and 76.1% of the total variance of RE among the four road types. The biplot of Figure 4 shows that ROR/S practically is the only variable that influenced dbRDA1 (given the high loading on this axis). The other two variables (R/US and R/BS) have a significant influence only on dbRDA2, which however explains only a very small share (0.7% of the fitted model and 0.5% of the total variation) of total variation of RE. Therefore, a reliable prediction model consisting of the ROR/S, R/US and R/BS variables may be adopted with satisfactory reliability, in order to predict, at least referring to the order of magnitude, the rill erosion rate in the analysed road types. This not only confirms the results of distance linear models applied in this study, but it is also in line with findings of the other authors; as a matter of fact [6,7,9,18], previously mentioned, have shown the importance of slopes of both road surface and upstream hillslope as well as bare soil presence as predictors of rill erosion in roads.

Conclusions
This study confirms the working hypothesis that the roads made of unpaved materials and, particularly, those in burned areas are more prone to rill erosion compared to rockpaved roads in areas subjected to wildfire. Moreover, a modeling approach to predict this erosion form could be based on geometric variables that are easy to be measured and processed, such as the slopes of both road surface and upstream hillslope as well as the percent bare soil over the road surface. However, the validity of these preliminary results must be confirmed in other experimental conditions, that is, in areas with different climate and soil conditions. Moreover, the modeling approach here suggested must be validated by larger applications of both field measurements and soil erosion prediction models. Other variables (e.g., rainfall erosivity, soil properties, traffic conditions, logging activities) must be considered to understand all drivers of soil erosion related to road design and construction and to identify the most effective methods for controlling and mitigating soil erosion rates. Moreover, the contribution of the most recent technologies (such as remote sensing techniques and image processing algorithms) to save time and money in mapping soil erosion sources in forest roads may be of help for future monitoring and management strategies. These methods are useful to support the activity of land managers in a better understanding of the magnitude of road erosion and developing efficient measures for its control and mitigation (such as road material or design used). The information obtained in this work could lead to a better design of forest roads that can effectively reduce rill erosion after wildfires.