Soil Trafficability Forecasting

This article introduces and evaluates a Soil Trafficability Model (STRAM) de-signed to estimate and forecast potential rutting depth on forest soils due to heavy machine traffic. This approach was developed within the wood-forwarding context of four harvest blocks in Northern and Central New Brunswick. Field measurements used for model calibration involved determining soil rut depths, volumetric moisture content, bulk density, soil resistance to cone penetration (referred to as cone index, or CI), and the dimensionless nominal soil cone index (NCI) defined by the ratio of CI over wheel foot print pressure. With STRAM, rut depth is inferred from: 1) machine dimensions pertaining to estimating foot print area and pressure; 2) pore-filled soil moisture content and related CI projections guided by year-round daily weather records using the Forest Hydrology Model (ForHyM); 3) accounting for within-block soil property variations using multiple and Random Forest regression techniques. Subsequent evaluations of projected soil moisture, CI and rut-depth values accounted for about 40 (multiple regression) and 80 (Random Forest) percent of the corresponding field measured values.


Introduction
With the wide-spread use of modern mechanized harvest machinery, forest and agricultural planners have to deal with the effects of heavy machine loads on soil rutting, compaction, erosion, and subsequent reductions in crop yields (Brady & Weil, 2008). Increased soil compaction reduces soil porosity, increases soil runoff, damages and crushes roots, and leads to reduced root growth due to decreased soil oxygen levels (Grigal, 2000;Horn et al., 2004;Bassett et al., 2005;Singer & Munns, 2006;Chen & Weil, 2011). Rutting results from combined soil compaction and soil displacement. Soil displacement occurs when soils within depressions and along slopes are moist to wet at or above field capacity (Raper, 2005;Naghdi et al., 2009). Depending on their slope orientations, ruts increase, decrease, or collect run-off (Sutherland, 2003;Antille & Godwin, 2013;Poltorak et al., 2018).
A measure of soil resistance to penetration, known as cone index (CI), is used to assess soil compaction and rutting. This measure varies from weak to strong as soil moisture content decreases, and this is particularly so in fine-textured soils. In contrast, sandy soils remain friable from wet to dry (Earl, 1997;Vaz, 2003;Dexter et al., 2007;Tekeste et al., 2008;Vaz et al., 2011;Kumar et al., 2012;Jones & Arp, 2017).
To comply with best forest management practices, various jurisdictions have established criteria as to what constitutes a rut. For example, the Province of British Columbia (Canada) classifies ruts > 15 cm deep to be a hazardous disturbance. The Pacific Northwest region of the USA defines ruts to be >15 cm deep as well (Page-Dumroese et al., 2000). The Province of Alberta considers any soil disturbance tracking deeper than 10 cm as ruts (Alberta Forest Products Association, 1994;Van Rees, 2002). A visual representation of minor to severe soil disturbance classes including rut depths is provided in Figure A1.
Earlier rut-depth prediction models by Maclaurin (1990) and Rantala (2001) related rut depth as caused by single wheel and single machine passes to the ratio of tire footprint pressure over soil resistance to penetration, also referred to as nominal cone index, or NCI (Table 1). Scholander (1974) and Saarilahti (2002a) extended this approach towards multiple passes. Vega-Nieva et al. (2009)   b u e = + + + + + D = wheel diameter; CI = cone index; n = number of passes; NCI = nominal cone index; CF = coarse fragment; M = pass-cumulative machine mass; a, a n , b 1 , b 2 , b 3 , p 1 , p 2 = multi-pass coefficients; r 1 = site specific calibration parameter; u m ; e m = random effects; NR: not reported. data for a sandy soil and a clay loam soil at varying moisture contents, and with a wood-forwarding rut-depth data by Meek (1996) suggestion to correct for coarse fragment content. A recent study by Sirén et al. (2019) related rut depth to number of machine passes, volumetric soil moisture content, and CI. All these studies, however, captured only a small portion of field-determined rut depth variations whether based on fixed or fixed plus random effects regression models. This article reports on the extent of improving wood-forwarding rut-depth modelling based on: 1) Using field-determined data for soil density, texture, organic matter, coarse fragment and weather-inferred ridge-to-valley soil moisture, CI, and machine-specific NCI variations for the rut-depth projection purpose.
2) Employing fixed and random forest regression techniques for optimizing these projections.
3) Selecting four of the eleven harvest blocks described in Jones and Arp (2019) for detailed analysis; this involved assessing wood-forwarding tracks across two shelterwood cutting blocks and across two commercial thinning blocks in Northern and Central New Brunswick, Canada.

Block Descriptions
Details regarding the wood-forwarding operations within the four harvest blocks chosen for this study are presented in Figure 1 & Figure 2 by location, and in Table 2 by block-specific attributes. These blocks are a subset of the 11 blocks used by Jones and Arp (2019) for their on-and off-track soil moisture and CI study.   Block 3 is a 20.9 ha shelterwood cut located in the southern reach of Notre Dame Mountains, within the Appalachians Mountain range. The area has a mean air temperature of 3.5˚C, with mean January and July temperatures at −12.9˚C to 17.6˚C, respectively. It has a mean precipitation of 1140 mm, with 310 mm as snow (Department of Environment and Climate Change Canada, 2016a). This block supports moderate to intolerant hardwood/mixedwood vegetation, mostly composed of red maple (Acer rubrum L.), yellow birch (Betula alleghaniensis Britt.) and balsam fir (Abies balsamea L.). Bedrock formations involve late Ordovician deep-water marine clastics, with topography ranging from gentle, rolling plateaus to steep valleys of hummocky ablation moraines.
Block 9 is a 25.5 ha shelterwood cut located in the Southern tip of the Miramichi highlands, roughly 50 km northwest from the Fredericton, the provincial capital for New Brunswick. The area has a mean air temperature of 5.5˚C, with mean January and July temperatures at −9.4˚C and 19.4˚C, respectively. Mean annual precipitation amounts to 1100 mm, with 250 mm as snow (Department of Environment and Climate Change Canada, 2016a). This block supports natural mixedwood and tolerant hardwoods, including Eastern hemlock (Tsuga canadensis L. Carrire), yellow birch, sugar maple (Acer Saccharum Marsh.), beech (Fagus grandifolia Ehrh.), balsam fir, Eastern white cedar (Thuja occidentalis L.), and black spruce (Picea mariana Mill). Surficial deposits vary from morainal sediments to ablation and boulder tills, underlain by early Devonian mafic volcanic rock and late Ordovician marine clastics.

Wood Forwarding
Three GPS-tracking wood-forwarding machines were used: a John Deere 1510E forwarder in Blocks 1 and 2, a John Deere 1110E forwarder in Block 3, and a Tigercat 635D in Block 9 ( Table 3). The resulting GPS data were processed to determine the number of passes per same track using point buffering and overlapping tools (Buja, 2012). Studies have shown that the number of wood-forwarding passes affects soil structure, compaction, and rut depth (Eliasson, 2005;Ampoorter et al., 2007;Farzaneh et al., 2012;Jones et al., 2018). Some studies have looked at mitigating multiple-pass soil rutting using brushmats (McDonald & Sexias, 1997;Labelle et al., 2015). Brushmats were not found along the forwarding trails in Blocks 3 and 9, but were present to varying degree in Blocks 1 and 2 as part of the commercial thinning operation. Studies have also shown that Open Journal of Forestry lower tire footprint areas due to increasing tire pressure increases rutting (Raper et al., 1995;Saarilahti & Antilla, 1999;Saarilahti, 2002b;Jun et al., 2004;Affleck, 2005;Sakai et al., 2008).

Data Sources
The data needed to process the temporal and spatial modelling of MC, CI, and soil rutting refer to: 1) Daily precipitation (snow and rain) and mean air temperature data, needed for the hydrological soil moisture calculations for Blocks 1-3 and Block 9. These data were obtained from the airport weather station records at Saint Leonard and Fredericton, respectively (Department of Environment and Climate Change Canada, 2016a).
2) Hydrometric stream discharge data, needed to calibrate ForHyM. These data were obtained for Blocks 1 -3 from the nearby Black Brook Watershed Research Site (2014), and for Block 4 from the Nashwaaksis stream data (Department of Environment and Climate Change Canada, 2016b).
3) LiDAR-generated bare-earth elevation data. These were downloaded from geoNB's website at 1 m resolution (GeoNB, 2015). 4) Digital soil maps (DSM layers) for soil texture (sand, silt and clay content), bulk density, coarse fragments, and soil organic matter for top 30 cm of soil, at 10 m resolution (Furze, 2019). 5) GPS-recorded machine-clearance pattern, spaced at 10 sec intervals along each wood-forwarding track. The generation of these data was described by Jones et al. (2018).

Field Measurements
Soil samples as well as volumetric soil moisture (MC v ) and cone index (CI) read-  Jones and Arp (2019). The samples were analyzed for texture (sand, silt and clay %), coarse fragments (%), organic matter content (%), and bulk density (g/cm 3 ). Wood-forwarding rut depths were determined at GPS-recorded intervals along the tracks. The number of passes along the tracks were obtained from the GPS-tracked wood-forwarding machine-clearance records (Jones et al., 2018) through counting the number of passes per same track using digital point buffering and overlapping tools (Buja, 2012).

Soil Moisture Modelling
The with log 10 DTW FIA referring to the logarithm of the rasterized depth-to-water index adjusted to the upslope flow initiation area (FIA) at the time of wood forwarding. In addition, OM DSM is the digitally derived soil organic matter raster (Furze, 2018; in %), DEM is elevation in meters, Track = 1 for ruts and 0 otherwise, and HW = 1 refers to locations dominated by hardwood, otherwise HW = 0. The values of model parameters required for Equations (1) and (2) are listed in Table 4.

Cone Index and Rut Depth Modelling
Rut locations and depths were GPS tracked and measured in the fall of 2014 where z is best-fitted field-measured rut depth, and where NCI adj = (NCI + 0.48 DTW FIA ) incorporates a further numerically derived block-by-block DTW adjustments.

Statistical Analyses
All analyses were performed in R (R Core Team, 2015) after combining the rasterized data layers into a covariate stack. A correlation matrix was created to show the general association pattern between the covariates, and the covariate association pattern so revealed was factor analysed. This was followed modelling MC PS , CI, and z n using multiple and Random Forest regression formulations (Breiman, 2001). The dataset for modelling MC PS (n = 394) and CI (n = 372) included all top 15 cm soil sampling locations inside and outside the tracks in Blocks 1 to 11. The dataset for modelling rut depth (n = 207) applied to Blocks 1, 2, 3, and 9 only.
The multivariate regression (MR) formulation via Equations (2) to (5)  MR and RF inputs comprised NCI (Equation (4)), CF DSM , OM DSM , depth to C horizon (C DSM ), and number of passes ("Passes"). The training process was set to generate and test 15 regression trees using 3 variables to split each tree node, followed by a 10-fold cross-validation process with 10 replicates (Kohavi, 1995).
The resulting model output produced field-determined versus Random-Forest fitted MC PS , CI, and z n scatter plots, and informed about the most MC PS -, CI-, and z n -predictive data layers. RF model performance was determined by evaluating the mean decreased accuracy and variable importance.
Both MR and RF techniques were used to determine which combination of variables would produced the best-fitting MC PS , CI, NCI and rut depth projections. With best-fitted results reported by listing the number of variables used, the root mean square error (RMSE), and the coefficient of determination R 2 .

Soil Trafficability Model (STRAM)
The information flow generated from modelling soil moisture, CI and rut depth was organized in the form of a Soil Trafficability Model framework (STRAM), as shown in Figure 3. Applying STRAM involved: 1) Initializing and calibrating ForHyM to determine MC PS for specific weather conditions, by block and stand type.

MCPS, CI, and Rut-Depth Analysis
The basic statistics of the plot-centered rut depth, MC PS , CF DSM , Sand DSM , OM DSM , D b,DSM , C DSM , elevation, DTW, CI, and NCI values are listed in Table 5. The correlation coefficients and their non-zero significance levels between these variations in association with block location are compiled in Table 6, along with the corresponding factor analysis. Factor 1 indicates that the rut-depth determinations related positively to number of wood-forwarding passes, to sand and organic matter content, and to soil depth (C DSM ), but negatively to NCI, as to be expected. In addition, rut depth was generally deeper in hardwood blocks (Blocks 3 and 9). Factor 2 reflects that the hardwood blocks at higher elevation and DTW locations had lower soil density and hence and lower soil resistance to penetration than the softwood blocks located at lower elevations.  The MR-generated results in Table 7 revealed that MC PS increased significantly with MC PS,DTW as influenced by wet weather and by decreasing DTW from ridge tops to low-lying areas next to water-filled flow channels. In addition, MC PS was higher inside than outside the tracks, and higher under hardwood vegetation but decreased towards the higher elevation blocks (Table 7). CI decreased with increasing pore-filled soil moisture content, but was noticeably higher within softwood blocks and tracks due to cumulative wheel-induced compaction. As per Equation (5), rut depth increased significantly with number of forwarding passes (log 10 Passes) and decreased significantly with increasing log 10 NCI. The non-linear DTW-adjustment for Equation (5) improved the MR regression results and related best-fitted scatterplots, but only to a small extent.
The RF results, also listed in Table 7, identified MC PS , DTW, elevation (DEM) and OM DSM as dominant MCPS influencing predictors. For CI, RF selected elevation, MC PS , Sand DSM , Depth, and track location served as dominant predictors. For rut depth, the RF process selected number of passes, NCI, DTW and CF as best predictors.
Compared to MR, the best-fitted RF-generated R 2 values for MC PS , CI and rut depth in Table 7 were considerably higher. The corresponding scatter plots in Figure 4 demonstrate the extent of the MR-versus RF-generated goodness-of-fit for each of these variables. In detail, the eight-times-out-of-ten conformance levels for MC PS , CI and rut depth respectively amounted to: ±14%, ±0.7 MPA and ±14 cm for MR, and ±4%, ±0.3 MPA, and ±4 cm for RF ( Figure 5).
The successive inclusion of additional predictor variables from most to the least significant was characterized by diminishing conformance gains, as shown in Figure 6. Using only one continuous predictor variable for MC PS led to an      variables. Having these variables as initial predictor variables was essential to address the otherwise unresolvable scatter that would otherwise be incurred by using continuous predictor variables only. Figure 7 shows the GPS-tracked wood forwarding tracks across Blocks 1, 2, 3, and 9 are overlaid on the delineated DTW < 1 m and hill-shaded DEM background, also shown are the number of passes (top) and field-measured rut depths (bottom) at each location (Figure 7). Among these blocks, Block 3 revealed a close association between deep rut depths and DTW, followed by Block 9 and Block 2.
Block 1 had no association between rut-depth and pass number. Rutting was deepest in Blocks 3 and 9 along multi-pass tracks across streams and wet-areas.   Open Journal of Forestry detailed in appearance and conformance than the MR projections, as to be expected from the scatter plots in Figure 4, and the results listed in Table 7.

MCPS, CI, and Rut-Depth Projections
Among the blocks, Block 9 proved to be the wettest, on account of 70 mm rain event two days before prior to the day of wood forwarding, and as reflected by the a consequential use of the DTW FIA = 0.25 ha data layer as dominant MC PS predictor via Equation (4). In contrast, Block 1 was found to be excessively dry, such that the DTW FIA = 16 ha projection was best to represent the field-determined MC PS values for this block outside the tracks. For Blocks 2 and 3, the MC PS was best presented using the DTW FIA = 1 ha assignment to reflect the May and June soil moisture conditions at the time of wood forwarding. Although both blocks had the same DTW FIA = 1 ha assignment, they differed in terms of CI-measured soil strength which was determined to be weaker for the hardwood block (Block 3) than for the softwood block (Block 2). Typically, shallow-rooting softwood forests grow on coarser and stonier soils with lower soil organic matter accumulations than deeper-rooted hardwood forests. Figure 10 presents an overlay of rut depth points on the corresponding MR (top) and RF (bottom) projections for Blocks 1, 2, 3, and 9 after 1, 10 and 50 passes, with better and more resolved RF than MR data-to-projection conformances. These plots confirm that number of passes and spatial variations in soil moisture are important rut-depth predictor variables. Actual rut depth, however, also depends on machine weight/load and soil physical properties, as quantified by way CI, NCI, and the variables listed in Table 7. Figure 10. Random Forest modelled 2-and 10-pass rut depths for Blocks 1, 2, 3, 9 with machines carrying full loads, with the field plot determinations overlaid, using FIA = 1, 16, 1, 0.25, respectively.

Discussion
While RF emulates field-measured values for MC PS , CI and rut depth considerably better than best-fitted MR values, it can only do so by systematically tracking those variables that best account for the overall data variations, including outliers. To be of general value, more testing is required to capture more of the MC PS , CI and rut depth variations across a wide range of soil and vegetation conditions. In this regard, the above regression results are at least consistent with general expectations. For example, moist to wet soils have low physical strength due to low particle cohesion (Kumar et al., 2012), and are therefore prone to traffic-induced compaction, displacement and rutting (Sutherland, 2003;Børgesen et al., 2006;Nikooy et al., 2016;Jones & Arp, 2017). To illustrate, Block 3 shows deeper ruts within the wetter DTW FIA = 1 ha marked area next to a stream. Block 1 with no significant rut-depth observations was cut following dry weather conditions during mid-June of 2014, with overall soil moisture levels best conforming to a DTW FIA = 16 ha flow-channel pattern. In contrast, deep ruts were encountered across Block 9 due to field operations in October 2014 following a 70 mm heavy rain event.
In terms of other physical soil properties, studies have shown that measured rut depth correlates positively with increased levels of OM in the soil (Sutherland, 2003;McFero Grace et al., 2006). In this study, rut depth also correlated positively with OM in Block 9 as revealed by the factor analysis results shown in Table 6.
Increased sand content generally contributes to low CI, NCI and therefore to increased rut depth, mainly due to low particle-to-particle cohesion (Balland et al., 2008;Brady & Weil, 2008;Kumar et al., 2012). In contrast, high coarse fragment content would increase CI, thereby decreasing rut depths. However, the overall CF variations within and across the blocks did not register this effect by way of MR, and only weakly so by RF. Typically, soils with high soil strength (high CI and NCI values) minimize soil disturbance (Antille & Godwin, 2013).
Significant and influential on the MR and RF results was the dependence of rut depth on the number of passes (p < 0.0001). In detail, this effect decreased non-linearly with increasing pass number due to gradually increasing soil compaction, as quantified above via Equation (5) (2018). This is especially so for high traffic areas such as wood landing sites, and along tracks involving a hundred passes or more (Carter et al., 2007;Taghavifar & Mardani, 2014;Jones et al., 2018).
Given the above moderate (MR) and strong (RF) data-to-projection conformances, it should be possible-at least in general-to project and forecast soil trafficability by following the flowchart given in Figure 3.  years. Generally, soil trafficability within this block would be best on solidly frozen ground, but worst during snowmelt periods when soils tend to be equally wet and partially thawed across the land. In detail, April and November would be the wettest months ( Figure 11). During spring, summer and fall, soil trafficability would vary by monthly variations in soil moisture, and by the reduction thereof through evapotranspiration, as primarily affected by precipitation, air temperatures and canopy leaf area. According to Figure 11, the driest summer For the spatial component of the soil trafficability projections, it is important to determine the likely upslope flow-channel initiation area for each block, as listed in Table 8. To some extent, these estimates would need to be modified by soil texture and coarse fragment content: FIA numbers should be higher for well drained and rocky soils, and lower for loamy and clayey soils. Extended droughts and frost periods would also increase FIA, therefore increasing the areas available for off-road traffic. However, even during winter, care needs to be given to not drive along or across flow channels with high upslope flow-accumulation areas. This care is needed because: 1) channels may not totally be frozen on account of upwelling seepage water, and 2) soil compaction and rutting along and across the channels would not only aggravate subsequent flooding but also soil and stream bank erosion. An example of year-round soil trafficability forecasting by month and related weather-imposed DTW FIA assignments is provided in Figure 12, with a focus on likely rut depths to be incurred by two wood-forwarding passes.  Perhaps the greatest impediment for correct soil trafficability forecasting is the lack of proper coarse fragment specifications, which should-ideally-represent total CF changes within and across soils. The CF data and data layers used for this purpose were, however, not revealed to be significant within-and across-block predictor variables for MC PS , CI and rut depth. In part, this was due to not including coarse fragments larger than gravel size in the soil sampling process.
However, where such sampling is sufficient, the following formulation for CI

Concluding Remarks
While this research demonstrates that the STRAM approach can be used to assess and project soil trafficability through coupling block-based soil surveys with temporal and spatial modelling techniques, there is a need for further data layer improvements using MR and RF modeling techniques. To this end, area-systematic rut surveys can now be conducted by, e.g., equipping off-road vehicles with GPS-tracking LiDAR-based rut depth sensors (Salmivaara et al., 2018). Similarly, Giannetti et al. (2017) proposed terrestrial portable laser scanners and Haas et al.
(2016), Pierzchała et al. (2016) and Launiainen et al. (2017) used unmanned airborne vehicles (UAVs) for rut-depth stereo imaging and evaluating underlying terrain conditions. In addition, advances in weather-affected digital soil property mapping will further assist soil trafficability mapping, and the validity of the same needs to be assessed systematically using area-wide post-operational rut-depth surveys. In this regard, MR and RF techniques will be helpful in terms determining how operationally induced soil compaction and rutting as observed and projected vary by topographic location, season and weather. Open Journal of Forestry

Appendix 2. Temporal Soil Moisture Modelling
The Forest Hydrology Model (ForHyM) uses daily temperature and precipitation (rain and snow) data as well as block-specific area, vegetation, soil horizon texture, depth, OM, and CF (Table A1) to predict soil moisture and temperature fluctuations through the soil (Arp & Yin, 1992;Yin & Arp, 1994;Jutras, 2012).
Calibrating ForHyM consists of comparing actual snowpack and hydraulic flow to modelled outputs and adjusting the output parameters (Table A2, Figure   A2).