ABSTRACT

This article introduces and evaluates a Soil Trafficability Model (STRAM) designed 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.

Keywords:

Soil Trafficability, Wood Forwarding, Plot Surveys, Regression Comparisons, Cartographic Depth-to-Water

1. 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) modified the NCI-based multi-pass regression model by examining the experimental

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₁, b₂, b₃, p₁, p₂ = multi-pass coefficients; r₁ = 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.

2. Study Locations, Materials and Methods

2.1. 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.

¹SWPL: softwood plantation, IntHW: intolerant hardwood, MW: mixedwood, TolHW: tolerant hardwood. ²BS: black spruce; WS: white spruce, BI: birch, YB: yellow birch; RM: red maple; SM: sugar maple; EH: eastern hemlock. ³CT: commercial thinning, CC: clear cut, SHW: Shelterwood, SC: select cut.

Blocks 1 and 2 are within 1 km of each other, and share similar site characteristics. They are located within the Western reach of the Chaleur uplands. Both blocks are white spruce (Picea glauca (Moench) Voss) plantations. Block 1 covers a two hectare area. Block 2 covers 20.8 ha, but was sampled acrossa 3.5 ha area. The region has a mean annual temperature of 3.8˚C, with mean January and July temperatures at −11.5˚C and 18.2˚C, respectively. Mean annual precipitation amounts to 1110 mm, with 300 mm as snow (Department of Environment and Climate Change Canada, 2016a) . Surficial deposits vary from morainal sediments to loamy lodgment tills, underlain by late Ordovician marine clastics.

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.

2.2. 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

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) .

2.3. 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) .

2.4. Field Measurements

Soil samples as well as volumetric soil moisture (MC_v) and cone index (CI) readings were obtained from the top 15 cm layer of undisturbed mineral soil across Blocks 1 to 11 after wood forwarding (Table 2) in the fall of 2014, as described in Jones and Arp (2019) . The samples were analyzed for texture (sand, silt and clay %), coarse fragments (%), organic matter content (%), and bulk density (g/cm³). 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) .

2.5. Soil Moisture Modelling

The Forest Hydrology Model ForHyM was used to estimate daily variations in pore-filled soil moisture content (MC_PS) on the well-drained soils in each block (Appendix). The MC_PS output so generated served to determine how MC_PS would have varied at the time of wood forwarding across the terrain of each block by way of season- and weather-adjusted cartographic depth-to-water (DTW) modelling (Murphy et al., 2009, 2011; White et al., 2012) .

The cross-terrain MC_PS projections were originally based on (Vega-Nieva et al., 2009; Jones & Arp, 2019) :

${\text{MC}}_{\text{PS,DTW}}=100-\left[100-{\text{MC}}_{\text{PS,top}}\right]\times {\left[\frac{1-\exp \left(-k_{mc}\text{DTW}\right)}{1-\exp \left(-k_{mc}{\text{DTW}}_{\text{top}}\right)}\right]}^{p_{mc}}$ (1)

where MC_PS,top refers to the pore-filled soil moisture content, DTW_top refers to at the highest elevation in each block, k_mc and p_mc are block-specific calibration parameters, with MC_PS,DTW = 100% along water-filled flow channels where DTW = 0 m. The MC_PS,DTW projections inside and outside wood-forwarding tracks were augmented via multiple regression analysis to account for changes in elevation, forest cover type (hardwoods, softwood, open areas), and soil organic matter content (OM_DSM), as follows:

$\begin{array}{l}{\text{MC}}_{\text{PS}}=70.88+0.26{\text{MC}}_{\text{PS,DTW}}-2.76{\log }_{10}{\text{DTW}}_{\text{FIA}}-0.98{\text{OM}}_{\text{DSM}}\\ -0.10\text{DEM}+10.22\text{HW}+4.09\text{Track};R^2=0.46\end{array}$ (2)

with log₁₀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.

2.6. Cone Index and Rut Depth Modelling

Rut locations and depths were GPS tracked and measured in the fall of 2014

immediately after wood forwarding, as listed in Table 2. Volumetric soil moisture, bulk density (D_b) and cone penetration measurements were obtained prior to these operations (Jones & Arp, 2019) . Once MC_PS was rasterized via Equations (2), CI (in MPa) and NCI were emulated using the following expressions (Jones & Arp, 2019) :

$\text{CI}=1.60-0.01{\text{MC}}_{\text{PS}}+0.03\text{Depth}+0.76\text{Track}+0.40\text{SW};R^2=0.41$ (3)

where Depth is soil depth (cm), Track = 1 signifies presence of track with Track = 0, is undisturbed soil, and SW = 1 signifies where softwoods are dominant, otherwise SW = 0. The normalized cone index NCI (dimensionless) was derived from:

$\text{NCI}=\frac{1000\text{CI}db}{W}\sqrt{\frac{\delta }{h}}\frac{1}{1+b/2d}$ (4)

with d as tire diameter (m), b as tire width (m), W as total wheel load (kN) given by (machine weight + load)/number of tires, h as section height of the tire (m), and $\delta =0.001\left(0.365+170/p\right)W$ as tire deflection, with p is tire inflation pressure (in kPa).

Relating rut depth to NCI and number of passes along the same track produced the following regression equation:

$z=\frac{40.8}{{\text{NCI}}_{\text{adj}}}\times {\text{Passes}}^{0.071\times {\text{NCI}}_{\text{adj}}};R^2=0.32$ (5)

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.

2.7. 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) served to explore the extent of functional dependence of MC_PS, CI, and z_n on weather, season, machine type, and digital soil mapping layers. The required input for modelling MC_PS via Equation (2) refers to MC_PS,DTW (Equations (1)), log₁₀DTW, OM_DSM, DEM, HW = 1 or 0, and Track = 0 or 1. The required input for modelling CI via Equation (3) refers to modelled MC_PS, Depth, SW = 1 or 0, and Track = 0 or 1. The required input for rut depth modelling via Equation (5) refers to modelled NCI and number of passes.

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.8. 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.

2) Projecting MC_PS across the terrain based on digital elevation, DTW, and soil mapping techniques, by block, using Equations (2).

3) Combining the weather-specific MC_PS projections with digitally generated soil properties to generate the CI, NCI and rut-depth data layers, by block, using Equations (3) to (5).

4) Optimizing these data layers through R-tool regression techniques (fixed effects, Random Forest) based on actual field observations where available.

3. Results

3.1. MC_PS, 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₁₀Passes) and decreased significantly with increasing log₁₀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² 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

RF-R² gain of about 0.80, while the corresponding MR-R² gain was limited to 0.2. For CI and rut depth, MR and FR had similar low initial R² values using, respectively, the binary 0, 1 Track and the grouped number of passes (“Passes”) 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.

Figure 8 provides a closer look regarding the extent of track rutting in Block 9 through the overlays of the DTW < 0.5 m delineated patterns on the hill shaded-DEM (top) and the surface-image (bottom) backgrounds. The surface image was generated a year after the wood-forwarding operations. At that time, rutting appearances had faded but remained prominent in the lower left corner of Block 9. Rutting > 40 cm deep occurred along multiple pass tracks where DTW < 0.5 m.

3.2. MC_PS, CI, and Rut-Depth Projections

Figure 9 presents the MR-and RF-generated projections and data points for MC_PS (top) and CI (bottom) for Blocks 1, 2, 3, and 9, with RF projections more

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. 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.

4. 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) and in conformance with Eliasson (2005) , Eliasson and Wästerlund (2007) , Botta et al. (2009) and Jones et al. (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. Fundamental to doing this is the combining of preceding and forecasted weather conditions (Appendix) with digital soil mapping. To this effect, Figure 11 provides a summary of daily summer through winter, weather conditions for Block 1 for a period of 14

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 occurred in 2005 followed by wettest fall. The wettest summers occurred in 2003 and 2011. Soil trafficability would also have been very poor in April 2005. In contrast, soil trafficability would have been best during the winters of 2004, 2011, 2013 and 2014 due to deep and prolonged snowpack and soil frost conditions.

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 may serve as general non-block specific CI predictor (Jones & Arp, 2017) :

${\log }_{10}\text{CI}=0.26-0.29\times \text{PS}-0.41\times {\text{MC}}_{\text{PS}}+1.04\times \text{CF};R^2=0.54$ (6)

where data for total coarse fragment content per soil volume are available, rut projections need to be forced towards zero when CF approaches 1. This can be done by, e.g. adjusting Equations (5) to

$z=\frac{40.8}{{\text{NCI}}_{\text{adj}}}\times {\text{Passes}}^{0.071\times {\text{NCI}}_{\text{adj}}}\times {\left(1-\text{CF}\right)}^2$ (7)

5. 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.

Acknowledgements

This research was supported financially through an NSERC Collaborative Research Project with J.D. Irving Limited as industrial partner. Special thanks go to Forest Watershed Research Centre at UNB, for facilitating this research, with additional support from Doug Hiltz for assistance with field sampling, and from Dr. Shane Furze for RF-analysis assistance and 10 m DSM layer provisions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Jones, M.-F., & Arp, P. (2019). Soil Trafficability Forecasting. Open Journal of Forestry, 9, 296-322. https://doi.org/10.4236/ojf.2019.94017

References

Appendix

Appendix 1. Rut Depth Severity Classes

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).