Locating and Projecting Hydrologically Induced Forest Road Vulnerabilities ()
1. Introduction
Constructing and maintaining forest roads in good condition (Swift, 1988; Ryan et al., 2004; Laschi et al., 2019) is a world-wide and costly challenge due to downhill water flow and subsequent roadbed erosion and sediment transfer from unpaved roadbeds (Keramati et al., 2020; Dodson, 2021; Lyons et al., 2023). Depending on traffic loads, slope, weather, season, and dry to moist and wet soil conditions, these roads deteriorate gradually through rutting, braiding, puddling and flooding. Through surveying, design, planning, maintenance, upgrading, and best management practices, these effects can be reduced (Aruga et al., 2007; Grace III & Clinton, 2007; Eastaugh et al., 2008; Zhao et al., 2022). Taking notes of these conditions ranges from on-the-ground observations and testing to remote sensing, including satellite imaging, UAV tracking, and terrestrial as well as airborne laser scanning (see, e.g. Akgul et al., 2017; Hrůza et al., 2018; Siafali & Tsioras, 2024). Studies about machine-effected tire loads onto roadbed penetrability have been reported by Suvinen et al. (2009), Vega-Nieva et al. (2009) and Jones and Arp (2017).
This article focuses on forest road trafficability and adjacent forest lands through digital elevation modelling (DEM) at 1-m resolution. This is done through: 1) image-based locating and indexing the type of road deterioration prior to road upgrading, and 2) applying this index through digital elevation modelling and mapping to the hydrological context of select road locations as affected by their topographic position and soil wetness in relation to slope, depression, and flow channeling. In so doing, this article builds on the DEM-based probability assessment of off-road machine-induced rut occurrences in clearcuts (Snow et al., 2024). Related trafficability studies have been reported by Niemi et al. (2017), Kankare et al. (2019), Salmivaara et al. (2020), Waga et al. (2020), Pundir and Garg (2022), Fjeld et al. (2024), and Kim et al. (2025). To this effect, Kankare et al. (2019) summarized “… that the risk of rutting damage is significantly affected by the local weather history, route planning, used machinery and the use of harvest residue to reinforce the tracks... but additional research is warranted to fully comprehend … the risk of rutting damage”.
2. Methods
2.1. Study Area
The case study refers to: 1) a 12 km section of the Aboujagane Road east of Moncton (Figure 1), and 2) is supplemented with a smaller 4 km section of the Penniac Road northeast of Fredericton, and an additional location in southwestern Nova Scotia. Currently, segments along these roads are mainly used for logging access and transport. For the Aboujagane Road section, elevations vary from 20 to 200 m north to south. Its area is located on red-grey sedimentary bedrocks of Pennsylvanian to Triassic age. This partially Ca-carbonate cemented formation is covered by shallow ablation and basal tills consisting of sandy to loamy texture with a low coarse fragment content. The soil, referred to as Parry soil (Colpitts, 1995), has a reddish color, compacts and ruts easily, and is prone to puddling and flooding along depressed and poorly drained road segments. For the Penniac Road section, elevations range from 70 to 180 south to north. Its area is covered by Harcourt soil that developed on a medium to coarse till on top of fine-textured basal till, all underlain by a red mudstone formation, also of Pennsylvanian to Triassic age (Colpitts, 1995), and compacts and ruts easily as well. The climate for both areas is maritime, with mean monthly temperatures ranging from −13˚C to 25˚C, and annual precipitation rates of about 1200 mm, with a quarter of that being snow. These areas were selected through image-inspecting the extent incurred road damage. The Aboujagane Road segment was used for model development. The Penniac Road and Nova Scotia locations served for image-based model verification.
![]()
Figure 1. Locator map for the Aboujagane and Penniac areas (top, red outlines), with close-ups (bottom) displaying the existing road networks. The depicted Aboujagane Road section (red line, bottom left) is used for model development. Areas A, B, and C are used for model verification along road segments, adjacent areas, and select power lines marked as red stippled lines in A and B. Background: GeoNB-ESRI Basemap.
2.2. Road Segment Inspections
Inspecting high-resolution LiDAR-generated DEMs generated in 2015 and 2017 and historical Google Earth and other imageries (GeoNB, ESRI World Images) from June 2015 to August 2024 revealed considerable variations in road segment quality along the Aboujagane Road and elsewhere. The best historical Google Earth Images (Figure 2) were used for marking and indexing 67 damaged locations (Figure 3) and 27 additional no-damage ridge-top locations along the Aboujagane Road. Doing so produced the road damage index (RDI) as an off-road trafficability indicator, symbolized as no damage 0, rutted 1, braided 2, puddled 3, and flooded 4, and also symbolized as 0, r, b, p, f, respectively, with examples in Figure 4. This evaluation was supplemented and confirmed: 1) by on-the-ground roadbed and culvert inspections and 2) AUV imaging in May 2020 (Figure 3). The resulting RDI values for the non-upgraded road segments were averaged for further analysis.
![]()
Figure 2. Best-focused historical Google Earth images centered on road flooding (Location 62) along the Aboujagane Road before (Sept. 2015, May 2017) and after (June 2022, Aug. 2024) road repair, centered on Location 62.
2.3. Terrain Analysis
The 2015 and 2017 LiDAR-generated Digital Elevation Model (NB DEMs; http://www.snb.ca/geonb1/e/dc/catalogue-E.asp) and the GeoNB data layers for forest and provincial roads, waterbodies, watercourses, and wetlands (http://www.snb.ca/geonb1/e/dc/catalogue-E.asp) were used to produce the following data layers via ArcGIS Pro processing:
1) Depressions and related depth (Ddepth, in m), determined by subtracting the original DEM from its depression-filled version.
2) Slope (in %), followed by focal 5-m radial smoothing.
3) Continuous flow direction (FD) and flow accumulation (FA, m2 converted to ha), using the D8 algorithm (Tarboton, 1997).
4) Log-transformed FA raster, i.e. logFAmax = log10(FAmax + 1), followed by focal 5-m radial smoothing around each raster cell to serve as a flow-accumulation index across each road point of interest.
Figure 3. Aboujagane road segment with 67 locations affected by rutting, braiding, puddling, and flooding. Left panel: Hill-shaded 1-m DEM, overlaid by DEM-derived flow channels with >4 ha upslope flow accumulation (white lines) and road-affecting watershed extent (blue overlay). Mid panel: GeoNB-registered wetlands and upgraded logging-road segments (red lines) with recent cut block activities as per August 2024 Google Earth imagery. Right panel: Close-up overlaid by the May 2020 UAV imaging track.
5) Flow channel lines with >0.1, >1, and >4 ha upslope FA.
6) Watershed borders for the FA > 4 ha flow lines where they cross the selected road segment.
7) Cartographically determined depth-to-water (DTW, in m; Murphy et al., 2011), calculated with slope as cost raster and with DTW = 0 where water can be expected to at the surface.
8) Topographic position index (TPI, in m; Weiss, 2001) based on the elevational difference between every DEM cell at center and its surrounding annulus 50 m away.
9) Probability of rut occurrences along single tracks, formulated as (Snow et al., 2024):
(1)
Figure 4. Road deterioration examples. Left panel: FA > 4 ha channel-induced flooding (f) at Location 39; FA > 0.1 ha channel-induced rutting (r) at Location 40. Right panel: FA > 0.1 ha channel-induced braiding (b) and puddling (p) at Locations 22 and 23, respectively. Thin and thick white lines: DEM-derived flow channels > 0.1 ha and > 4 ha.
The workflow associated with establishing these layers is presented in Figure 5. The hill-shaded version of the DEM and the FA > 0.1 ha raster layers were subsequently used: 1) to locate the hydrological centers of each of the 94 RDI locations, and 2) to extract the values for Ddepth, FA, 5-m mean focal logFAmax, mean focal 5-m Slope, DTW, TPI, and PRrut at the RDI points.
2.4. Statistical Analysis
The 94 data points so generated were subject to multivariate regression analysis, with the averaged RDI numbers as dependent variable and Elevation, Slope, Ddepth, DTW, logFAmax, TPI, and PRrut as independent variables. This was done stepwise to determine which of the resulting models and associated RDI predictor variables would contribute the most to the road-incurred damage rating, from no rutting to flooding. Once established, the best-fitted results were further refined by identifying and removing model-associated RDI outliers.
2.5. Mapping and Model Verification
The best-fitted model so generated was used to map RDI across the Aboujagane and Penniac study areas in Figure 1. The verification involved inspecting the semitransparent RDI overlays on: 1) the hillshaded DEMs for Areas A, B, and C, and 2) the corresponding Google Earth images, all overlaid by the DEM-generated flow channels with FA > 4, >1, and >0.1 ha. Image-revealed damage locations along the road segments and powerlines so selected were circled to visually ascertain their correspondence to green (0), yellow (2) to red (4) coloration of the RDI map. This verification was further extended by checking the RDI map performance for a forest road location in southwestern Nova Scotia (Latitude 44.271; Longitude 64.872).
Figure 5. Workflow for the data layers used for assessing terrain conditions along road segments and adjoining areas of interest.
3. Results and Discussion
3.1. Road-Damaging Attributes
The image-determined RDI numbers and their averages along the non-repaired road segments are listed in Table 1 by location. image year, and associated FA, Slope, and Ddepth numbers. As shown in Figure 6 for May 2017, puddling and flooding were most noticeable within depressed water-accumulating road locations, while rutting and braiding occurred along lower slopes. The corresponding post-repair RDI evaluation based on the August 2024 imagery (Locations 26 to 60) is also presented in Figure 6, with RDI values mostly rated zero, except for rutting and puddling at Locations 45 and 46 next to the re-aligned roadbed. At Location 51, post-repair rutting also occurred but only as per the June 2022 image.
3.2. RDI Analysis Results and Interpretations
The stepwise regression analysis selected Ddepth, Slope, logFAmax, and PRrut as the
Table 1. Road damage index (RDI) numbers along the pre-repaired Aboujagane road segments grouped by upslope FA classes: >4 ha, 1 to 4 ha, 0.1 to 4 ha, and <0.1 ha. Indexed as no damage 0, rutted 1, braided 2, pooled 3, and flooded 4 by image year, averaged across years, and with FA, slope, and depression depth also noted for each location where applicable.
Figure 6. Elevation profile of the damaged road locations (Locations 1 to 65), color-coded no damage (0), rutted (r), braided (b), pooled (p), and flooded (f), based on the May 2017 and June 2022 Google Earth images. Black arrows point to RDI = 0 locations.
main RDI predictors. The correlation among these variables varied from −0.83 (Ddepth vs. Slope) to 0.633 (logFAmax vs. PRrut). The remaining coefficients varied from 0.158 to 0.374. The best-fitted results so obtained are listed in Table 2 from Models 1 to 6. This stepwise progression is based on and affected by the following considerations:
1) Log-transforming FA via logFAmax linearizes the FA effect along and across roads.
2) Using the 5-m focal mean Slope serves to reduce the slope variations across the DEM-profiled roads and adjacent ditches and bulldozed materials where present.
3) The DEM-assessed depression depths depend on the extent of water filling at the time of the LiDAR-based DEM scanning process, being naturally deepest when dry and equal to 0 when flooded.
4) Cells with low (0) to high (1) PRrut values are easily spotted by overlaying the model-generated PRrut raster on the high-resolution images and DEMs for the area of interest.
Table 2. Best-fitted multivariate regression results generated through stepwise reduction of least significant DRI predictor variables; n = 94 (Model 1), 93 (Models 2, 3, and 4), 88 (Models 5 and 6).
Model |
RDI Predictors |
Regression Coefficients Summary |
Coeff. |
Std. Error |
t-value |
p-value |
R2 |
RMS |
1 |
Intercept |
1.6 |
±0.24 |
6.5 |
<0.0001 |
0.819 |
0.625 |
PRrut |
1.8 |
±0.22 |
9.0 |
<0.0001 |
Ddepth |
4.62 |
±0.92 |
5.0 |
<0.0001 |
Slope |
−0.169 |
±0.031 |
−5.5 |
<0.0001 |
logFAmax |
0.365 |
±0.068 |
5.3 |
<0.0001 |
2 |
Intercept |
0.83 |
±0.22 |
3.8 |
0.0003 |
0.771 |
0.703 |
PRrut |
2.70 |
±0.21 |
12.9 |
<0.0001 |
|
Ddepth |
5.23 |
±1.02 |
4.3 |
<0.0001 |
|
Slope |
−0.146 |
±0.034 |
−4.3 |
<0.0001 |
3 |
Intercept |
0.10 |
±0.15 |
0.7 |
0.515 |
0.729 |
0.761 |
PRrut |
2.48 |
±0.22 |
11.3 |
<.0001 |
Ddepth |
6.0 |
±1.1 |
5.4 |
<.0001 |
4 |
Intercept |
0.98 |
±0.25 |
3.9 |
0.0001 |
0.704 |
0.795 |
PRrut |
3.13 |
±0.22 |
14.4 |
<0.0001 |
Slope |
−0.176 |
±0.038 |
−4.6 |
<0.0001 |
5 |
Intercept |
0.07 |
±0.15 |
0.1 |
0.644 |
0.720 |
0.76 |
PRrut |
3.15 |
±0.21 |
15.0 |
<0.0001 |
6 |
PRrut |
3.23 |
±0.12 |
27.6 |
<0.0001 |
0.719 |
0.772 |
The results in Table 2 were generated after excluding outlier locations as follows: from Model 1 (Loc. 66), from Models 2, 3, 4, and 5 (Loc. 66, 67), and from Models 5 and 6 (Loc. 12, 20, 21, 49, 66).
Progressing from Model 1 to Model 6 slightly decreased the R2 values and slightly increased the RMSE values while the RDI coefficients remained near the best-fitted Model 1 values. Removing logFAmax caused the least changes. Removing Slope and then also Ddepth rendered the best-fitting RDI scattergram (Figure 7) less uniform because PRrut was designed to only map probable single-track off-road rut occurrences, regardless of whether rutted, braided, puddled, or flooded. Nevertheless, PRrut still emerged as the strongest RDI predictor because removing PRrut from Model 1 caused the largest best-fitted R² drop, i.e. from 0.819 to 0.632, with no outliers removed. The wider area-influence of PRrut compared to the road-focused Ddepth and Slope determinations is likely due to the fact that PRrut captures the soil-rutting effect in water-accumulating depressions and along permanent to ephemeral flow channels quite well, as quantified via Equation (1) and as shown in Snow et al. (2024).
![]()
Figure 7. Image-determined versus best-fitted RDI scattergrams for Model 1 after removing its Location 66 outlier (red dot, left) and Model 6 (right) after removing its outliers at Location 12, 20, 21, 49, 66, and 67, i.e. dots marked red, yellow and green in the left panel.
The reason that logFAmax is least significant among the RDI predictors relates to the fact that rutting, braiding, puddling, and/or sustained flooding is much reduced at culvert-protected road locations for which FA tends to be >4 ha, as was the case at Locations 20, 41, 48, 52, and 66. In contrast, 1) culvert washout had occurred at Locations 11 and 27 where FA = 65.5 and 76.7 ha, respectively, and 2) roadbeds were fully wetland-flooded at Locations 39 and 62.
The installed culverts were all corrugated and had a diameter of D = 500 mm. According to the Manning Pipe Flow formula (https://en.wikipedia.org/wiki/Manning_formula), this would accommodate even-flow stream discharge rates up to:
(2)
and assuming that culvert roughness is represented by η = 0.022, and culvert slope (Sculvert) is set at 0.5% to allow for fish-passage. For even-flow run-off discharge resulting from a ppt = 100 mm/day storm when soils are fully saturated, this rate would be exceeded at all road-stream crossings where
(3)
Such occurrences would lead: 1) to cumulative roadbed soil loss, as noted for Locations 48 and 52, where FA = 11.4 and 13.9 ha, and/or 2) to washouts as noted for Locations 18 and 39 where FA = 161 ha and 104 ha, respectively. Noting that the discharge-accommodating culvert radius can be estimated from:
(4)
It follows that Locations 18 and 39 would require culvert diameters amounting to 1260 and 1070 mm to prevent overflow and potential washout in response to an even-flow 100 mm/day storm event. However, single-storm discharge peaks at about three times higher than the corresponding even-flow discharge expectations (Gray, 1962). Hence, Locations 18 and 39 would need culvert diameters amounting to 1900 and 1600 mm, respectively, to avoid 100 mm/day caused overflow and washout occurrences. Similarly, the corresponding diameters for Locations 48 and 52 would amount to 700 and 750 mm, respectively.
3.3. RDI Mapping and Verifying
Using Model 1 in Table 2 for the RDI mapping purpose matches the image-observed and circled road-damage locations quite well. In Figure 8, this is illustrated for the 38 to 43 Aboujagane Road Locations. At these locations, the May 2017 image-observed RDI variations revealed:
1) Channel-flow induced braiding (38, RDI = 2.8).
2) Flooding (39, RDI = 3.4).
3) Rutting on flat but slightly depressed uphill terrain (40, RDI = 1.1).
4) Rutting and puddling due to roadbed flow channeling (41, RDI = 3.4; 42, RDI = 3.1; 43, RDI = 3.3).
5) Wetland-expanded flooding at Location 39.
6) No damage (0) along nearby uphill locations.
By expressing Model 1 in equation form:
(5)
It can be noted that RDI: 1) can take on values larger than 4, especially so with increasing values for Ddepth and logFAmax, and 2) can become negative with increasing Slope values. If so, RDI > 4 would indicate a yet greater possibility of road flooding. In contrast, RDI < 0 would imply no damage, as RDI projected along an upgraded and non-upgraded Aboujagane Road segment in Figure 9.
Figure 8. Comparing road-located rutting, braiding puddling, and flooding along the Aboujagane Road Locations 38 to 43. Top and bottom: Google Earth images from May 2017 (prior to repair) and August 2024 (after road repair). Middle: Hill-shaded RDI traffic colors indexed green to red from no damage (0) to flooded (4).
Extending the Model-1 RDI projections, region-wide revealed a general correspondence between image-observable and model-projected patterns regarding traffic-induced road and trail damage. Figures 10-14 demonstrate this not only for Sections A, B, and C in Figure 1, but also for a similarly distant location in
Figure 9. DEM-derived hill-shaded RDI raster (Model 1, top) and May 2017 Google Earth image (bottom) focused on active and inactive logging roads at and near Locations 22 to 25, with DEM-derived flow channels (thick lines FA > 4 ha; thin lines FA> 0.1 ha) overlaid.
southwestern Nova Scotia. In these illustrations, road damage generally occurs within the red RDI-indicated damage coloration along and/or next to adjacent non-active logging road segments, and also along powerlines. In contrast, this coloration tends to remain green where roadbeds were restored prior to LiDAR scanning, except where the red road coloration prevails at or near FA > 4 ha road-stream crossings. This suggests potential vulnerability to road rutting, braiding, puddling, flooding, and washout at these locations.
Reviewing Figures 8-14 indicates that the mapped RDI patterns can also be used to: 1) optimize road locations away from traffic-compromising slopes and water-accumulating areas to assess off-road rutting potentials, and 2) plan on- and off-road operations as affected by weather-, soil- and slope-anticipated conditions to avoid rutting. In this regard, rutting and flooding can be expected to occur in the lower-lying and water-accumulating sections of the red areas. As a further guide, it is recommended to overlay the DEM-projected RDI maps with officially delineated stream, wetland and waterbody maps. The advantages of
doing so include: 1) ensuring that already mapped waterbodies and wetlands become RDI > 4 entities, and 2) systematically checking and correcting for potential flow-blocking road-stream crossings.
It should be noted that the RDI mapping process as described above does not inform when newly constructed or upgraded forest roads succumb to rutting, braiding, puddling, and flooding. It does, however, inform where these roads are likely to become vulnerable due to factors such as:
1) Road washout at road-stream crossings with large upslope water-accumulating areas.
2) Culverts that are either too small, partially to completely plugged, or have incurred shape deformations (Mai et al., 2014).
3) Water-logging and flooding along roads next to poorly draining depressions.
Figure 10. Circled image-observed road and powerline damage within Section A (Figure 1) west of the Aboujagane Road Location 67, overlayed on: 1) the May 2017 Google Earth image (top) and 2) the hill-shaded Model-1 RDI raster (bottom). Also shown are the DEM (2017) generated FA > 0.1 and > 4 ha flow channels, represented by thin and thick white lines, respectively.
Figure 11. Circled image-observed road and powerline damage locations west of the Aboujagane Road in Section B of Figure 1, overlayed on: 1) the May 2017 Google Earth image (top) and 2) the hillshaded Model-1 RDI raster (bottom). Also shown are the DEM-generated FA > 0.1, > 1, and > 4 ha flow channels as represented by thin white lines, blue lines, and thick white lines, respectively.
Figure 12. Circled image-observed road and powerline damage locations along the Penniac logging road segments in Figure 1 (Area C), overlayed on: 1) the April 2021 Google Earth image (top), and 2) the hillshaded Model-1 RDI raster (bottom). Also shown on bottom, the DEM-generated FA > 0.1, > 1, and 4 ha flow channels represented by thin white lines, blue lines, and thick white lines, respectively. White circles refer to the RDI and FA > 0.1 ha identified road damages. Latitude 46.010; Longitude 66.497. Yellow rectangle: Figure 14.
Figure 13. Close-up of Figure 13, marking puddles (white circles) and eroded or partially eroded roadbed locations above road-stream crossings with FA > 1 ha (red rectangles) on the April 2021 Google Earth image (top) and hillshaded Model-1 RDI raster. Also shown FA > 0.1, > 1, and > 4 ha DEM-generated flow channels (bottom). Latitude 44.271; Longitude 64.872.
Figure 14. DEM-generated RDI and rut probability pattern for a location in Nova Scotia (middle), showing a close match between 1) the June 2019 Google Earth image (bottom), 2) the logged area, and 3) the associated road conditions. Logging continued until July 2021, leading to deeper rutting along the harvest tracks near the DEM-generated FA > 1 ha flow channels and adjacent RDI pattern (top panel). Latitude 44.271, Longitude 64.874; approximately 200 km south of Area A in Figure 1.
4. Concluding Remarks
Due to the DEM-generated RDI predictor variables, it is now possible to: 1) determine how forest road conditions vary over time, and 2) verify this where multi-year surface imagery of sufficient quality is available, road segment by road segment. Doing so by way of overlays reveals: 1) how road segments would be affected by moist to wet roadbed conditions, and 2) where roads would be vulnerable to flooding and washouts. Where surface images of high resolution and with non-obscuring clouds and/or tree shadows and snow cover are not available, on-the-ground inspections and UAV-based imaging would be needed for map-verifying road damage indexing. In this regard, the above DEM-based Aboujagane RDI model also proves to be applicable elsewhere with similar hydro-topographic forest conditions, as demonstrated by the examples in Figures 12-14. This is in part implied by using the DEM-generated PRrut variable that, by itself, was found useful for DEM-mapping and image-checking potential rutting risks within cut blocks and also along trails and powerlines (Snow et al., 2024).
Acknowledgements
Financial support from a 2021 NSERC CRD grant and related J.D. Irving Ltd. project collaboration regarding “Comprehensive soil mapping at 1-m resolution with operational forest management applications” is gratefully acknowledged.