Estimating soil moisture and subsequent resistance to cone penetration is at the base of forecasting potential soil disturbance effects due to off-road machine traffic. For example, modern forest harvesting operations including wood forwarding can lead to substantial soil compaction, rutting and displacements, rut-induced water logging and re-direction of flow patterns leading to soil erosion, operation inefficiencies, and increased wear of machinery component, especially when these operations are not properly timed. Machine-induced soil compaction and associated rut-induced soil displacements commonly occur on moist to wet ground, are long-lasting, and affect the growth of remaining or planted vegetation (Cambi et al., 2015; Solgi et al., 2018). Impacts are strongest along trails with multiple wood-forwarding passes, and on wood landing sites (Jones et al., 2018). To remain productive, post-harvest soils need to remain well drained with soil bulk densities at <1.5 g∙cm−3 (Soane & van Ouwerkerk, 1994; Sutherland, 2003; Bassett et al., 2005; Brady & Weil, 2008; Chen & Weil, 2011).
Forecasting soil compaction and penetrability across time and terrain is, however, difficult due to weather-affected soil moisture conditions, and meter-by-meter changes in soil substrates and properties (Elbanna & Witney, 1987; Smith et al., 1997; Vaz et al., 2001) pertaining to soil texture, coarse fragment content, and organic matter content, with penetrability-affecting soil cohesion increasing with decreasing particle size. The presence of soil organic matter modifies this effect due to organically-supported soil aggregation and related pore-space stabilization. The presence of coarse fragments reduces soil compaction and penetrability by increasing the force needed to displace these fragments downward and laterally (Baetens et al., 2009; Rücknagel et al., 2013). Soil penetrability leads to deep rutting on wet soils mainly due to soil displacement and on moist soils mainly due to soil compaction. Dry and dense soils are mostly resistant to penetration, compaction and rutting. Across terrains, soil penetrability changes due to:
1) changes in soil moisture content, which generally varies from well to excessively well drained on ridge tops to moist to water saturated soils along streams, shorelines, wetlands, in low-lying depressions, and in toe-slope seepage zones;
2) changes in soil type, depth, texture, coarse fragment and organic matter content, and bedrock exposure;
3) changes in tracked versus non-tracked ground, which in turn depends on extent of machine foot print, load, number of passes, and surface conditions as affected by the presence of stumps, roots, logging slash, rocks, depth of forest litter, and brushmats.
This article focuses on presenting a framework to model temporal and spatial weather-, forest-, terrain-, and machine-induced plot-by-plot changes in volumetric soil moisture (MCV) and soil penetrability cone index (CI) for a case study in northwestern New Brunswick, Canada. This study involved transect-sampling eleven forest blocks that were subject to clear cuts, selection cuts, shelterwood cuts, and pre-commercial thinning at different times of the year. The plot-generated MCV data were expressed in terms of pore-space filled soil-moisture content (MCPS) to serve as a useful CI predictor (Vega-Nieva et al., 2009). The data so obtained were used to model MCV, MCPS and CI spatially and temporally. The spatial modelling component addressed emulating the data variations across the terrain from ridge tops to valleys as proposed by Vega-Nieva et al. (2009). The temporal component involved emulating the data variations by weather and season using the Forest Hydrology Model (ForHyM, Jones & Arp, 2017).
The goodness-of-fit of the resulting best-fitted MC and CI models was evaluated in terms of uncertainty, and non-randomness indicators. This was done in two ways: 1) focusing on the field-determined data only, and 2) determining the extent to which province-wide elevation, cartographic depth-to-water (DTW), and soil property data layers (Furze, 2018) could be used for block-specific MC and CI projections, and hence, soil trafficability prediction purposes. Knowing where and when locations are subject to rutting has become a vital component of planning best forest management practices. To that effect, rut-avoidance regulations are already in place across many federal and provincial jurisdictions. For example, ruts > 15 cm deep are deemed to represent hazardous soil disturbances (Alberta Forest Products Association, 1994; Page-Dumroese et al., 2000; Van Rees, 2002).
2.1. Block Descriptions
The harvest blocks for this study were selected across the mid-western to northwestern sections of New Brunswick (Figure 1). Table1 informs about block-specific attributes pertaining to species composition, harvesting type with and without brushmats laid down along tracks, soil type, elevation, slope and aspect. Information about soil association, landform, and lithology for each block can be found in TableA1 (Appendix).
Northwestern Uplands (NWU)
Blocks 1 - 5 are located within the Southern Uplands Ecoregion of New Brunswick. Mean annual air temperature is 3.6˚C. Mean monthly temperatures vary from −5.3˚C (January) to 12.5˚C (July). The area has a mean precipitation of 1140 mm with 310 mm as snow, inclusively (Department of Environment and Climate Change Canada, 2016a). The forested vegetation mainly consists of sugar maple (Acer Saccharum Marsh.), balsam fir (Abies balsamea L.), yellow birch (Betula alleghaniensis Britt.), and black spruce (Picea mariana Mill.). Blocks 1, 2, 4, 5 involved white and black spruce plantations, whereas Block 3 involved naturally grown birch and maple trees. The bedrock formations of the area consist of late Ordovician deep water marine clastics. The topography comprises gently rolling till-covered plateaus with steeply incised valleys.
Midwestern Uplands (MWU)
Blocks 8, 10 and 11 are located on the Miramichi Caledonia Highlands, which are also part of the Southern Uplands Ecoregion of New Brunswick. Mean annual air temperature is 5.2˚C. Mean monthly temperatures vary from −3.5˚C (January) to 13.9˚C (July). The area has mean precipitation of 1180 mm with 280
Figure 1. Digital elevation model for New Brunswick (10 m resolution) with regional upland and lowland delineations, block locations, and weather station locations.
Table 1. Block description by geographical location, forest type, stand properties (elevation, slope, aspect), and soil association type, as further described in Table A1.
1SWPL: softwood plantation, IntHW: intolerant hardwood, MW: mixedwood, TolHW: tolerant hardwood; 2BS: black spruce; WS: white spruce, BI: birch, YB: yellow birch; RM: red maple; SM: sugar maple; EH: eastern hemlock; 3CT: commercial thinning, CC: clear cut, SHW: Shelterwood, SC: select cut. 4Tracks brushmatted.
mm as snow (Department of Environment and Climate Change Canada, 2016a). Block 8 consisted of a naturally regenerated tolerant hardwood stand, with mostly sugar maples and yellow birch trees. The bedrock formations of the area consist of early Devonian felsic plutons, generally overlain by loamy lodgment tills and sandy glaciofluvial outwash sediments.
Blocks 6, 7 and 9 are located in the midwestern portion of the Continent Lowlands Ecoregion of New Brunswick. Mean annual air temperature is 5.5˚C. Mean monthly temperatures vary from −2.8˚C (January) to 13.8˚C (July). Annual precipitation amounts to 1100 mm, with 250 mm as snow (Department of Environment and Climate Change Canada, 2016a). Blocks 6, 7 and 9 supported natural mixedwood and tolerant hardwood stands, comprised of Eastern hemlock (Tsuga canadensis L. Carrière), yellow birch, sugar maples, beech (Fagus grandifolia Ehrh), and some balsam fir, Eastern white cedar (Thuja occidentalis L.), and black spruce. The bedrock formations of the area consist of early Devonian mafic volcanic and late Ordovician deep-water marine-clastics, generally overlain by ablation and boulder tills, glaciofluvial sediments and eskers.
2.2. Data Sources and Processing
Data layers needed for the spatial and temporal evaluation and modelling the plot-by-plot data involved:
1) producing region-wide digital elevation models (DEMs, 1 m resolution) from GeoNB’s LiDAR elevation point cloud data using LAS tools;
2) producing the cartographic depth-to-water layer (DTW), for the purpose of emulating plot-by-plot changes on soil moisture content as affected by weather conditions at the time of field sampling (see below);
3) using the province-wide 10 m resolution soil property data layers for soil bulk density (Db), coarse fragments (CF), organic matter (OM), and texture developed by Furze (2018); this development used topographic, climatic and geological data layers to emulate soil drainage, horizon depth, depth, texture, coarse fragment content, organic matter content, and bulk density data as specified for 12,058 geo-referenced soil pedon locations;
4) obtaining daily weather records for daily rain and snow amounts and snowpack depth from weather station, within or near the NWU, MWU and LL regions (Department of Environment and Climate Change Canada, 2016a) for block-specific soil moisture emulations;
5) obtaining daily stream discharge records from hydrometric monitoring stations representative of stream water flow within or near the NWU, MWU and LL regions for hydrological model calibration (Department of Environment and Climate Change Canada, 2016b) to ensure that the soil moisture emulations were consistent with regional weather and stream discharge events;
6) acquiring forest inventory and road layers, to navigate to and access harvest blocks.
All raster and shapefile data layers were assembled with ArcMap software, using the same projection system (NAD 1983 CSRS New Brunswick Stereographic). The transect plots within blocks were georeferenced. The resulting coordinates were used to extract plot-specific data values from the elevation, DTW, Db, CF, OM, and soil texture data layers. The DEM layer was also used to determine mean elevation, slope, and aspect for each block (Table 1), needed as block-specific soil-moisture modelling input.
2.3. Field Measurements
Soil property evaluations were done along 696 geo-referenced transect plots inside and outside wood-forwarding tracks 1) for Blocks 1 to 9 intermittently from May 27, 2014 through November 21, 2014, and 2) for Blocks 10 and 11 in June 2013 (Table 1). The plots within blocks were established pairwise, each pair 100 m apart containing one plot within and one plot adjacent to the track on undisturbed soil, 10 m apart (Figure 2).
Soil samples were retrieved from the top 15 cm of mineral soil. Each sample was placed into labeled freezer bags for storage. Samples were dried in a forced-air oven at 75˚C for 24 hours, then crushed and passed through a 2 mm sieve to separate the fine earth from the CF. The latter was used to determine CF%. The former was used to determine 1) sand, silt and clay for each sample (Sand% + Silt% + Clay% = 100%), using the hydrometer method (Shelrick & Wang, 1993),
Figure 2. Penetrometer being used inside and outside forwarder tracks without brushmats.
and 2) soil carbon% by oven-dry weight using a LECO CNS-2000 analyzer. The resulting soil carbon numbers were converted into soil organic matter via Equation (1).
OM = 1.72 × Carbon(%) (1)
Soil density samples were also collected from each plot by tapping and extracting a metal cup of known volume (85 cm3) into the mineral soil. Following frozen storage, these samples were weighed and dried in a forced-air oven at 75˚C for 24 hours. From this, the non-sieved soil Db also containing roots and coarse fragments was determined as per Equation (2).
The fine-earth bulk density of the soil between rocks and roots was estimated using the oven-dry weight fractions of sand (SandW), organic matter (OMW) and mineral soil depth (cm) as predictor variables as per Equation (3) (Balland et al., 2008).
A Humboldt digital cone penetrometer (cone base = 1.5 cm2; cone angle 60˚) was used to determine soil penetrability through recording CI in MPa at 15, 30, 45, and 60 cm depths within each sampling plot. Similarly, a Delta T HH2 moisture meter (TDR: a time-domain reflectometer) was used to determine the volumetric soil moisture content of the fine-earth fraction between coarse fragments for each plot at 15 cm mineral soil depth. Five CI and MCV readings were obtained for each plot: one at the center point (location of GPS coordinates), and four arranged at each cardinal direction (north, south, east, and west), 1 m apart from center. For blocks 10 and 11, the CI measurements were taken once per subplot to the deepest depth reached with the penetrometer. Measurements were not recorded if obstructed by surface rocks, roots or logs. The subplots were then averaged to provide single numbers per plot for further analysis.
Since CI tends to be related to pore-filled soil moisture content (MCPS, see Vega-Nieva et al., 2009), it was important to infer soil pore space (PS) and MCPS from the soil particle and bulk densities as follows (Balland et al., 2008):
with organic and mineral particle densities set to DOM = 1.3 and DMin = 2.6 g∙cm−3, respectively.
2.4. Temporal MC Variations
Soil moisture was modelled block-by-block over time using the temporal aspatial Forest Hydrology Model (ForHyM) (Arp & Yin, 1992; Yin & Arp, 1994; Jutras, 2012). This model uses 1) daily weather records for temperature, and precipitation for input (Table 2), 2) block-specific specifications for topography (elevation, slope, and aspect), 3) soil-surveyed properties (horizon depth, texture, depth, OM, CF) as per mapped soil type (Table 3), and 4) vegetation type and % canopy closure. These specifications were needed to emulate daily soil moisture, temperature, snowpack conditions and stream discharge. In this model, soil permeability at saturation was empirically related to Dp, Db, and Sandw by setting
; R2 = 0.80 (7)
In ForHyM, Ksat and snowpack are further calibration-adjusted to quantify the extent of inflow (lateral flow) along layers versus downward percolation (vertical flow) into layers (Table 2, FigureA1). In soils where soil density increases with depth, downward Ksat is general less than the Equation (7) specifications. Also of note is the ForHyM formulation for water retention at field saturation (FCW), given by
Table 2. ForHyM calibration variables for snowpack and saturated soil permeability by block, including weather and hydrometric station locations for primary ForHyM input data.
1 Department of Environment and Climate Change Canada (2016a); 2Black Brook Watershed Research site (2014).
Table 3. Soil profile information used to initialize ForHyM each block.
S = sand, SL = sandy loam, SCL = sandy clay loam, L = loam.
; R2 = 0.96 (8)
where FCW, SPW, and SandW all refer to total weight fractions per dried and gently crushed fine-earth soil that passes through a 2 mm sieve. In combination, any ForHyM generated block-by-block soil moisture output by soil layer is, apart from daily weather records for rain, snow and air temperature and topography specifications is also a function of and soil bulk density and sand, organic matter and coarse fragment content. Specifically, Equation (8) implies that soil moisture retention at field capacity increases with increasing organic matter content and decreasing bulk density, sand and coarse fragment content (Balland et al., 2008).
2.5. MC and CI Variations, Plot-by-Plot
In general, soil moisture increases from ridges and steep slopes to water-saturated wetlands, depressions and hyporheic zones adjacent to streams, rivers, lakes and shores. To some extent, this tendency is quantitatively related to the cartographic depth-to-water index (DTW). This index represents the least elevation rise to any point on the land away from the nearest open water locations, where DTW = 0. Conceptually, soil moisture levels should therefore decrease as DTW increases. For the changes in volumetric and pore-filled soil moisture content over time, this can be expressed as follows (Vega-Nieva et al., 2009):
with MCPS, DWT and MCV, DWT as pore-filled and volumetric soil moisture contents at any point in the landscape in relation to DTW, with MCPS0 as pore-filled moisture content on top of the ridge, where DTW = DTWmax, and with DTW = 0 at open-water locations such as stream and lakes. The local flow-channel network was generated by applying the D8 flow-accumulation algorithm to the 1 m resolution LiDAR-generated digital elevation model for New Brunswick, as described in Figure 3. The seasonality extent of the resulting flow network was emulated by varying the minimum upslope flow-accumulation area for flow initiation (termed flow initiation area, or FIA for short) by season and weather, as illustrated in Figure 4.
In Equation (9), parameters kmc and pmc (0.5 and 1.5, respectively) quantify how MCPS,DTW varies with DTW across areas of interest as MCPS0 varies from, for example, 20% on ridges to 100% at and near water-filled flow channels. When soils become uniformly saturated, MCPS,DTW = 100% everywhere. At other times MCPS0 ≤ MCPS,DTW < 1.
Figure 4. Cartographic DTW ≤ 1 m pattern by minimum upslope open-channel flow initiation area (FIA), i.e., DTWFIA, to emulate changes in soil moisture content by season and weather. For most of New Brunswick and the forested areas across Canada, the DTW ≤ 1 m pattern for FIA = 4 ha reflects end of the summer soil moisture and drainage conditions from very poor to moderate.
Since soil moisture at any point in the landscape are also affected by block-by-block and plot-by-plot variations in soil properties and wood-forwarding tracks, MCPS also needed to be evaluated as function forest cover and soil properties inside and outside wood-forwarding tracks. This was done by setting:
with Forest Cover (SW, HW, and MW), Block, Track and Brushmat coded 1 when applicable and 0 when not. The elevation variable accounted for additional topographic variations as they existed block-to-block across the NWU, MWU and LL regions. The soil properties refer to plot-by plot variations in soil texture (sand, silt, clay content), Db, CF, OM, and layer depth. The block-by-block and plot-by-plot CI variations were evaluated similarly by setting:
2.6. Statistical Analysis, Followed by Generalized Block-Generated MCPS and CI Projections
The block- and plot-descriptive, field, laboratory, ForHyM-generated and raster-extracted data were compiled into a single datasheet, were summarized, and were subjected to correlation, factor, multivariate regression analyses in R (R Core Team, 2015), using MCV, MCPS and CI as dependent variables. The resulting regression models for Equations (9)-(11) were assessed in terms of scatterplots for the actual versus modelled values, the best-fitted regression coefficients for each independent variable of non-zero significance values (p-values < 0.01), and R2 and RMSE goodness-of-fit indicators. In addition, the conformance levels between actual and model projected 0 - 10, 10 - 20, 20 - 30, … % MCPS classes and 0 - 0.5, 0.5 - 1, 1 - 1.5,… MPa CI classes were evaluated in terms of confusion matrices, non-randomness, and cumulative conformance probabilities to determine the uncertainty range of best-fitted model projections. Subsequently, the MCV, MCPS and CI regression evaluations were repeated using variables amenable for area-wide projection purposes. These variables referred to 1) DTW, elevation, digitally-generated province-wide soil property rasters for Sand, CF, OM and Db (Figure 5; Furze, 2018), 2) block-, and season- or weather-specific assignments for FIA, 3) soil depth, and 4) forest cover type for non-tracked soil conditions. The process of doing so by soil depth and season (month) is illustrated in Figure 6.
3. Results and Discussion
3.1. General Observations and Trends
The plot-based determinations for sand, silt, and clay, OM, CF, MCPS, and CI are summarized per block in Figure 7 by box plots and in terms of number of plots, and minimum, maximum, average and standard deviation values in TableA2 (Appendix). These plots reveal that the MCPS, CI and Db were not always higher inside than outside the wood-forwarding tracks, as one would expect (Allen, 1997; Han et al., 2009; Labelle & Jaeger, 2011; Cambi et al., 2015; Solgi et al., 2018). Possible reasons for the lack of systematic MCPS, Db, and CI increases refer to 1) the wide range of operational multitrack observations across ground conditions comprised of rocks, stumps, and variable forest litter depths including brushmats, and 2) differences in MCPS, Db, and CI surveying methodologies, e.g., using Db and MC sensors that required no soil displacement (Labelle & Jaeger, 2011) and hydraulic cone penetrators (Cambi et al., 2015; Solgi et al., 2018). The results reported below were obtained through manual soil extraction and probe insertions.
In terms of tracks with and without brushmats, measured values for Db, CI, and MCPS for the top mineral soil 15 cm were somewhat lower and less variable
Figure 5. Spatial raster examples for DTW (1), DEM (2), Db (3), CF (4) for block-specific MCPS projections (Block 3).
Figure 6. Digital soil mapping process combining topographic and soil features  with temporal features (B) to generate best-fitted weather-affected MCPS  and CI  projections for each block (e.g., Block 3) by way of best-fitted Equations (9)-(11) models.
Figure 7. Boxplots CI, MCPS, CF, Sand and Db within the top 15 cm soil layer inside and outside the wood-forwarding tracks, by Block.
under matted than non-matted tracks (Figure 8). Systematic increases for Db and CI from the matted and non-matted tracks as reported by Han et al. (2009) and Labelle et al. (2015) were therefore not found, likely due to strongly varying soil, root and rock conditions along the tracks.
Figure 8. Scatterplots pertaining to CI, Db, and MCPS inside and outside matted and non-matted tracks.
In terms of the higher trending MC levels on matted versus non-matted plots, Roberts et al. (2005) and Moroni et al. (2009) reported lower cumulative soil moisture losses under matted than non-matted tracks likely through mulch-related shading, insulation, and lowered evaporation benefits.
The overall trends between MCPS and CI versus CF, sand, OM and Db are presented in Figure 9. As shown, increasing CF generally lowered MCPS presumably due to 1) increasing porosities between the fragments, 2) subsequent increases in soil moisture loss due to increased soil permeability, and 3) enhanced evaporation due to fact that CF conduct heat faster than soil (Poesen & Lavee, 1994; Chow et al., 2007).
The increasing MCPS trend with increasing Db relates to decreasing pore space. Similarly, decreasing MCPS with increasing OM% relates to increased pore space due to OM-facilitated soil granulation (Han et al., 2006). There is also the possibility of decreasing MCPS due to increasing hydrophobicity as OM-containing soils dry out (Vogelmann et al., 2013). In contrast, the increasing MCPS trend with increasing Sand% may be due to preferential pore space saturation in sandier topsoil portions.
In terms of CI, increasing levels of CF contribute to increasing soil resistance to penetration because of skeletal soil stabilization (Manuwa, 2012), and the need to push coarse fragments to the side to gain greater penetration depths. Increasing Sand% generally decreases soil strength,, due to lower particle-to-particle cohesion. The opposite occurs with increasing Clay% (Raven et al., 1999). The trend between increasing CI and increasing soil depth (Figure 10) is due to depth-related increases in Db (Equation (3); Carter et al., 2007).
Analyzing the trends among the MCPS and CI affecting variables more closely produced the correlation matrix and its factor analysis results in Table 4. The bolded Factor 1 - 3 loadings indicate that:
1) Factor 1 is positively associated with increasing CI, CF, Db, elevation and tracks, but negatively associated with increasing soil organic matter content.
2) Factor 2 is positively associated with higher elevations where tolerant
Figure 9. Scatterplots for top 15 cm MCPS (top) and CI (bottom) versus top 15 cm CF, sand, OM and Db.
Figure 10. Boxplots of field-determined CI versus soil depth for Blocks 1 to 9, including significance value of implied trend. Blocks 10 and 11 do not have results by depth. Also shown (bottom): summary boxplots of all CI determinations off-track and inside tracks.
Table 4. Pearson’s correlation matrix (below 1,1,1... diagonal) with significance levels (above 1,1,1... diagonal) for each variable pair. Also shown: factor analysis results following oblique rotation.
hardwood forests dominate, where brushmats were rarely used, and where soils are somewhat sandier.
3) Factor 3 is positively associated with increasing MCPS but negatively associated with increasing CF, Sand, and decreasing DTW, as to be expected. Also, soil pores tend to be drier at higher elevations but fill up more easily and remain wet or moist longer inside than outside the tracks.
3.2. Temporal Soil Moisture Derivations
The ForHyM-generated results for daily soil moisture, stream discharge, snowpack depth and depth of soil frost are illustrated in Figure11, which also displays daily weather records (rain, snow snowpack depth, air temperature) for Blocks 1 and 9 from 2011 to end of 2014. The daily MCPS0 output so produced block-by-block (Table 5) was obtained through calibrating the ForHyM output with reported snowpack depth and cumulative stream discharge values for each of the three regions (FigureA1). These calibrations involved adjusting the ForHyM parameters for snowpack depth and lateral versus vertical layer-by-layer soil permeability, as detailed in Table 5.
Fitting the actual plot-specific MCPS determinations for the top 15 cm of mineral soil versus the corresponding log10DTW values with FIA = 1 ha across all the plots by way of Equation (9) produced the following equation once FIA and MCPS0 was specified for each block:
This result improved considerably by using log10DTW1ha in combination with MCPS,DTW:
In detail, the inclusion of log10DTW spread the MCPS prediction range from about 30% to 75% [Equation (12)] to about 20% to 90% [Equation (13)]. Even further improvements were obtained by including the plot determinations for Db, Elevation, HW and Track as independent variables. This led to the best-fitted MCPS regression results listed in Table 6, and to Equation (14):
Figure 11. Daily variations in Block 1 (left) and Block 9 (right) for air temperature and precipitation (ForHyM input; top) and stream discharge, MCPS for top 15 cm of mineral soil, estimated field capacity, snowpack depth, and frost depth (ForHyM output; bottom). Regarding ForHym details, see Jones et al. (2018).
Table 5. Mean block-specific soil property values including soil moisture content at ridge top (MCPS0) and minimum upslope open-channel flow initiation area FIA.
The corresponding best-fitted MCV results, also listed in Table 6, led to:
Table 6. Best-fitted MCPS% and MCV% models based on using plot-generated [Equation (14) and (15); top] versus plot-projected [Equation (16) and Equation (17); bottom], listing significant regression variables and their coefficients, standard error estimates, and t- and p-values, together with R2, RMSE values and sample size (n).
Units: MCPS, DTW and OMDSM in %; DTW and Elevation in m; Db in g/cm3; HW Blocks and Tracks: 1 when present, and 0 otherwise.
Based on the t-value entries in Table 6, one determines that the significance contributions of all the MCPS and MCV predictor variables follow these sequences:
for MCPS: Db ≈ MCPS,DWT ≈ Elevation ≈ HW > log10DTW > Tracks > Sand > CF;
for MCV: MCPS,DWT ≈ Elevation ≈ HW > log10DTW > Sand > Tracks > Db ≈ CF.
As to be noted, there is a high to low Db regression coefficient and significance change from MCPS [Equation (14)] to MCV [Equation (15)] due to the greater dependence of MCPS on Db via Equation (5) and Equation (6) while the significance levels for the predictor variables remain about the same. Also note that the best-fitted R2 and RMSE values drop from MCPS to MV due to the narrower MCV data range.
The + and − signs for the regression coefficients in Equation (14) and Equation (15) follow the expected trends for MC, namely:
1) MC decreases with increasing DTW but increases as MCPS0 (and therefore MCPS,DWT) increase, as to be expected.
2) MC decreases towards higher and generally steeper elevations.
3) MC is higher inside than outside soil tracks.
4) Soils under hardwood cover tend have higher MC values than elsewhere, in part due to deeper rooting and higher top-soil OM accumulation.
5) The contributions of Sand and CF (both negative) to MC were also found to be significant, but only marginally so. Excluding these from the analysis changed R2 for MCV from 0.51 to 0.50, and for MCPS from 0.63 to 0.60 (details not shown).
The best-fitted MCPS and MCV equations―obtained from substituting the plot-generated values for MCPS,DWT with the block-assessed FIA values and with plot-level Equation (3) estimates for OM generated from the digital soil layers for OM% (labelled OMDSM)―are as follows (Table 6):
The corresponding actual versus best-fitted scatterplots associated with MCPS via Equations (12)-(14) and Equation (16) are shown in Figure 12.
Using Equation (14) produced the block-by- block presentation in Figure 13, which also shows the field-determined MCPS plot values for the top 15 cm of soil outside and inside the tracks. Visually, there is a general agreement between the off-track field determinations and the corresponding projections, with off-track MCPS generally higher than inside track.
Grouping the off- and in-track projections and corresponding field determinations into 10% MCPS classes produced the MCPS confusion matrix and associated MCPS cumulative conformance plots in Figure 14. In summary, the ±5% MCPS conformance level so found amounts to 34% across the 20 < MCPS < 100% range, but increases to 81% and 96% as the conformance uncertainty increases to ΔMCPS ≤ ±15 and ±25%, respectively. The possibility of randomly achieving
Figure 12. Scatterplots of actual versus best-fitted MCPS regression results. Using Equation (12) and (13) demonstrates the extent of block-specific DTW projections on MCPS. Using Equation (14) demonstrates the improvements obtained by adding field-determined soil density, elevation, forest cover type, and track versus non-tracked specifications to the MCPS regression analysis. Using Equation (16) demonstrates that reasonable MCPS can also be generated from digitally mapped variables (DTW, OMDSM) and block-specific assignments.
these conformance levels have kappa values of 0.26, 0.73 and 0.92, respectively. The kappa value for random chance agreements equals zero, by definition.
3.3. CI Derivation
The best-fitted CI results using plot-determined and plot-projected data are listed in Table 7 and led to the following equations:
where 0.03 log10DTW − 0.02 MCPS0 in Equation (19) replaced −0.01 MCPS in Equation (18). These results affirm that the plot-determined CI values, as to be expected, increase with decreasing soil depth, are higher inside that than outside
Figure 13. Overlay of field- determined MCPS values inside and outside tracks on the Equation (16) generated MCPS projections for Blocks 1 to 11. Note: Blocks 9 - 11 have outside-track MCPSvalues only.
Figure 14. Confusion matrix for actual versus Equation (13) projected MCPS (left) and best-fitted cumulative conformance probabilities for Equations (12)-(14), and Equation (16) (right). Blue squares align 1-on-1 on actual versus projected MC% classes, while the shades from dark red to white repersent the number of 1-on-1 class-by-class projection differences.
the tracks, are higher under softwood forest cover, decrease with increasing soil moisture content, and increase with increasing soil density. The associated significance levels vary as follows:
for Equation (18): Track > Depth > MCPS > SW > Db
for Equation (19): Track > Depth > SW ≈ DTW1ha > MCPS0
Table 7. Best-fitted CI models based on using plot-generated [Equation (18)] versus plot-projected [Equation (19)] and log10 plot-generated CI [Equation (20)], listing significant regression variables and their coefficients, standard error estimates, and t- and p-values, together with R2, RMSE values and sample size (n).
Other variables such as Sand, Clay, CF and OM would also affect CI, but these variables were not found to make additionally significant CI contributions to the plot-determined or projected values. Quantitatively, Equation (18) indicates that a change in MCPS from 0% to 100% would lead to a change in CI by 3.1 MPa (all other conditions remaining the same). In comparison, CI would on average increase by 2.6 MPa from 0 to 1m soil depth. In comparison, CI values are about 0.7 MPa higher inside than outside tracks.
The scatterplots associated with Equation (18) and Equation (19) are presented in Figure 15, with the former more heteroscedastic than the latter. Log-transforming CI the best-fitted results as listed in Table 7 produced an even scatterplot in Figure 15, based on the following equation:
The overlays of the field-determined on the Equation (18) projected CI values for the top 15 cm of soil are shown by Block in Figure 16, which also shows the corresponding values at 30, 45, and 60 cm depth for Block 3.
Grouping the field-determined and projected CI values into 0.25 MPa classes produced the CI confusion matrix and associated cumulative conformance plots in Figure 17. As shown, the CI confusion matrix has an actual with projected
Figure 15. Scatterplots of actual versus best-fitted CI scatterplots associated with Equations (18)-(20). Left: Using Equation (18) demonstrates the additional improvements obtained by adding forest cover type, soil depth and track versus non tracked specifications to the CI regression analysis. Middle: Using Equation (19) demonstrates that CI estimates can be generated from block-specific and weather-dependent MCPS0 and DTW determinations. Right: Using Equation (20) demonstrates how the actual log10CI values relate to the corresponding log10CI projections.
Figure 16. Equation (18) CI projections at 15 cm depth for Blocks 1 - 9 outside wood-forwarding (top), and for 15, 30, 45 and 60 cm soil depths outside (left) and inside (right) wood-forwarding tracks in Block 3 (bottom).
Figure 17. Confusion matrix for modelled vs Equation (18) projected CI (left) and best-fitted cumulative conformance probabilities for Equation (18) and Equation (19) (right). Blue squares align 1-on-1 on the actual versus projected MC% classes, while the shades from dark red to white represent the number of 1-on-1 class-by-class projection differences.
ΔCI ≤ 0.25 MPa class conformance of 32% across the 0.5 < CI < 4.5 MPa range. This conformance increases to 78% and 95% as the uncertainty range for CI increases from ±0.5 and ±0.75 MPa, respectively. Achieving these agreements has kappa values equal to 0.23, 0.67 and 0.88, respectively. In comparison with the MCPS kappa values, this suggests that the CI projections are slightly more random than the MCPS projections.
3.4. Forecasting Block-Specific Soil Moisture and Penetrability
The above results show that the combination of emulating layered soil properties and daily soil moisture simulations produces non-random soil moisture and cone penetrability projections across harvest blocks of varying topography, forest cover, substrates and weather conditions. To this effect, the extent of pore-filled moisture content can be predicted with a tolerance of ±15%, 8 times out of 10 within the 20% to 100% soil moisture range. Similarly, the cone penetrability index can be predicted with a tolerance of ±0.5 MPa, 8 times out of 10 within the 0.2 to 4.5 MPa range.
In summary, generalizing similar results for the purpose of soil-trafficability forecasting block-by-block requires:
1) a 10 m resolution DEM layer with a tolerance of at least ≤ 2 m, 8 times of out 10 (Furze et al., 2017);
2) a digital soil mapping process that emulates the required soil property layers at a general conformance level of 80% (Furze, 2018);
3) an assessment of the season and weather dependent minimum upslope open-channel flow initiation area (FIA) as this could vary from ≤0.25 to ≥8 ha (Figure 4);
4) a forest hydrology model that can be used to estimate weather- and season-affected soil moisture content (MCPS0) for ridge-top soils, to initiate, e.g., the Equation (15) and Equation (18) calculations (Jones & Arp, 2017), based on block-by-block elevation, slope, aspect, texture, Db, OM, CF, vegetation type and % canopy closure.
The MC and CI results so produce are projected in Figure 18 for hardwood
Figure 18. ForHyM modelled MCPS for ice and water combined, in relation to soil moisture retention at field capacity (top) with spatially projected MCPS [Equation (16)] and CI [Equation (19)] maps on block 3 (hardwood forest and soil conditions; middle) and block 4 (softwood forest and soil conditions; bottom).
and softwood areas somewhat larger than Block 3 and Block 4. This was done for weather conditions that varied from essentially frozen in February to wet in May, dry in July, and moist in October. The blue shading in these plots suggests where rutting would occur on account of elevated MC and lowered CI values once wood-forwarding operations cross the areas so marked. For the most part, areas prone delineations as these vary by weather and season. Least rutting would have occurred in February (ground mostly frozen) and in July and August (ground mostly dry), but would have been most extensive in May and late fall (ground mostly moist to wet).
3.5. Additional Comments
While the above study provides a useful framework for delineating and evaluating soil trafficability restrictions due to changes in soil moisture and soil penetrability, there is a need to validate the generality of this framework by:
1) testing across all dominant and minor soil types within the same region;
2) repeating the same across other regions;
3) conducting sequential MC and CI determinations as ground conditions transit from wet to dry and from frozen to non-frozen;
4) advancing the digital soil property mapping process from coarse textured DEMs (Furze, 2018) to 1 m LiDAR-generated terrain elevation resolution. Doing so will likely increase the significance by which projected sand, clay, organic matter, coarse fragments and bulk density gain significance as MC and CI predictors.
Repeated topsoil moisture and density testing should ideally be done without probe insertions, i.e. via nuclear moisture-density gauge (Corns, 1988; Brown et al., 1998; Labelle & Jaeger, 2011). CI testing would further benefit from hydraulic rather than manual testing (Hooks & Jansen, 1984; Herrick & Jones, 2002). Doing so would assist in reducing manually generated errors, inconsistencies, and biases, but would increase sampling expense.
There are other DEM-utilizing soil MC projection techniques. Most notably among these is the topographic wetness index (TWI; Beven & Kirkby, 1979; Sørensen et al., 2006). This index relates soil wetness to local upslope flow-accumulation areas and slopes, i.e., TWI = ln (flow accumulation area/slope). Comparative studies revealed that TWI-inferred MC depends strongly on TWI-cell averaging (Murphy et al., 2009a, 2011; Ågren et al., 2014), while DTW-inferred MC can directly be related to field-based DTW measurements (Dobbie & Smith, 2006; Murphy et al., 2009b).
Optimal TWI-based MC interpretations are generated by filling noise-generated DEM pits, followed by changing the focal search radius for cell-centered mean elevation differences, slope and TWI (Southee et al., 2012). The resulting MC correlations with TWI had R2 values ranging from 0.06 to 0.36, with best-results obtained for the 5 to 10 m cell-size range. Attempts focused on analyzing soil moisture conditions from spectral surface images tend to produce mixed results, especially for soils covered by vegetation and vegetation litter (Njoku & Entekhabi, 1996; Das & Paul, 2015). Best results are generally observed for bare to open mineral soil surfaces (Jackson, 1993; Engman & Chauhan, 1995; Wang & Qu, 2009; Das & Paul, 2015).
Further digital soil property mapping improvements can likely be obtained using LiDAR-generated DEMs at 1 m resolution. Doing so would likely increase the significance by which projected sand, clay, organic matter, coarse fragments and bulk density gain significance as MC and CI predictors.
With reasonable soil- and weather-informed MC and CI projections, it is feasible to generate machine- and load-specific rut depth maps as affected by machine type, loads, tire dimensions, and number of wood-forwarding passes, as demonstrated by Vega-Nieva et al. (2009) and Jones et al. (2018). In turn, such projections can be used to evaluate existing traffic-induced soil disturbance impacts as demonstrated by Campbell et al. (2013).
The above soil MC and CI assessment is limited to forested areas within northwestern New Brunswick, Canada. The approach taken emulates block-specific soil MC and CI variations with R2 = 0.61 and 0.40, respectively. The dominant predictor variable refers to daily soil moisture levels as affected by weather, DTW, sand, Db, CF, soil depth, stand type (HW or SW) and sampling location, inside or away from harvest tracks. Other potentially important MC and/or CI contributing variables yet to be addressed refer to the presence and extent of stumps and roots, depth of local forest litter accumulations, local variations in micro-topography pertaining to mounds and pits, and the number of passes per track. Implicitly addressed are the variations in soil OM due to the associations between OM, Db and soil pore space via Equations (3) to (6). The effect of brush-matting on MC and CI was found to be too variable to make significant contributions to the best-fitted MC and CI regression results. The lack thereof is likely related to branch-specific soil compaction and displacement at the centimeter scale.
This research was supported by an NSERC-CRD supported Forest Trafficability project, with further support from J.D. Irving, Limited. Special thanks go to the Forest Watershed Research Center at UNB for field sampling: Doug Hiltz, Mark Castonguay, Clara Ndekuringe Dennis, and Tanner Sagouspe.
Figure A1. Modelled vs measured ForHyM generated calibration outputs for snowpack and stream discharge for Northwestern Uplands (top), Midwestern Uplands (middle) and Lowlands (bottom).
Table A1. Soil association, landform, and lithology for Blocks 1 - 11.
Table A2. Soil physical and moisture properties by block.
1Blocks 10 and 11 measured max CI per plot.
 Ågren, A. M., Lidberg, W., Strömgren, M., Ogilvie, J., & Arp, P. A. (2014). Evaluating Digital Terrain Indices for Soil Wetness Mapping—A Swedish Case Study. Hydrology and Earth System Sciences, 18, 3623-3634. https://doi.org/10.5194/hess-18-3623-2014
 Allen, M. M. (1997). Soil Compaction and Disturbance Following a Thinning of Second-Growth Douglas-Fir with a Cut-to-Length and Skyline System in the Oregon Cascades. MSc Thesis. Corvallis: Oregon State University.
 Arp, P. A., & Yin, X. (1992). Predicting Water Fluxes through Forest from Monthly Precipitation and Mean Monthly Air Temperature Record. Canadian Journal of Forest Research, 22, 864-877. https://doi.org/10.1139/x92-116
 Baetens, J. M., Verbist, K., Cornells, W. M., Gabriels, D., & Soto, G. (2009). On the Influence of Coarse Fragments on Soil Water Retention. Water Resources Research, 45, 1-14. https://doi.org/10.1029/2008WR007402
 Balland, V., Pollacco, J. A. P., & Arp, P. A. (2008). Modeling Soil Hydraulic Properties for a Wide Range of Soil Conditions. Ecological Modelling, 219, 300-316.
 Bassett, I. E., Simcock, R. C., & Mitchell, N. D. (2005). Consequences of Soil Compaction for Seedling Establishment: Implications for Natural Regeneration and Restoration. Austral Ecology, 30, 827-833. https://doi.org/10.1111/j.1442-9993.2005.01525.x
 Brown, D. S. P.-D. R. E., Jurgensen, M. F., & Mroz, G. D. (1998). Comparison of Methods for Determining Bulk Densities of Rocky Forest Soils. Soil Science Society of America Journal, 63, 379-383.
 Cambi, M., Certini, G., Neri, F., & Marchi, E. (2015). The Impact of Heavy Traffic on Forest Soils: A Review. Forest Ecology and Management, 338, 124-138.
 Campbell, D. M. H., White, B., & Arp, P. A. (2013). Modeling and Mapping Soil Resistance to Penetration and Rutting Using LiDAR-Derived Digital Elevation Data. Journal of Soil and Water Conservation, 68, 460-473. https://doi.org/10.2489/jswc.68.6.460
 Carter, E. A., Aust, W. M., & Burger, J. A. (2007). Soil Strength Response of Select Soil Disturbance Classes on a Wet Pine Flat in South Carolina. Forest Ecology and Management, 247, 131-139. https://doi.org/10.1016/j.foreco.2007.04.026
 Chow, T. L., Rees, H. W., Monteith, J. O., Toner, P., & Lavoie, J. (2007). Effects of Coarse Fragment Content on Soil Physical Properties, Soil Erosion and Potato Production. Canadian Journal of Soil Science, 87, 565-577. https://doi.org/10.4141/CJSS07006
 Corns, I. G. W. (1988). Compaction by Forestry Equipment and Effects on Coniferous Seedling Growth on Four Soils in the Alberta Foothills. Canadian Journal of Forest Research, 18, 75-84. https://doi.org/10.1139/x88-012
 Dobbie, K. E., & Smith, K. A. (2006). The Effect of Water Table Depth on Emissions of N2O from a Grassland Soil. Soil Use and Management, 22, 22-28.
 Elbanna, E. B., & Witney, B. D. (1987). Cone Penetration Resistance Equation as a Function of the Clay Ratio, Soil Moisture Content and Specific Weight. Journal of Terramechanics, 24, 41-56. https://doi.org/10.1016/0022-4898(87)90058-9
 Furze, S., Ogilvie, J., & Arp, P. A. (2017). Fusing Digital Elevation Models to Improve Hydrological Interpretations. Journal of Geographic Information System, 9, 558-575.
 Han, H.-S., Page-Dumroese, D., Han, S.-K., & Tirocke, J. (2006). Effects of Slash, Machine Passes, and Soil Moisture on Penetration Resistance in a Cut-to-Length Harvesting. International Journal of Forest Engineering, 17, 11-24.
 Han, S.-K., Han, H.-S., Page-Dumroese, D. S., & Johnson, L. R. (2009). Soil Compaction Associated with Cut-to-Length and Whole-Tree Harvesting of a Coniferous Forest. Canadian Journal of Forest Research, 39, 976-989. https://doi.org/10.1139/X09-027
 Herrick, J. E., & Jones, T. L. (2002). A Dynamic Cone Penetrometer for Measuring Soil Penetration Resistance. Soil Science Society of America Journal, 66, 1320-1324.
 Hooks, C. L., & Jansen, I. J. (1984). Recording Cone Penetrometer Developed in Reclamation Research. Soil Science Society of America Journal, 50, 10-12.
 Jones, M.-F., & Arp, P. A. (2017). Relating the Cone Penetration and Rutting Resistance of Soils to Variations in Soil Properties and Daily Moisture Variations over Time. Open Journal of Soil Science, 7, 149-171. https://doi.org/10.4236/ojss.2017.77012
 Jones, M.-F., Castonguay, M., Jaeger, D., & Arp, P. A. (2018). Track-Monitoring and Analyzing Machine Clearances during Wood Forwarding. Open Journal of Forestry, 8, 297-327. https://doi.org/10.4236/ojf.2018.83020
 Labelle, E. R., & Jaeger, D. (2011). Soil Compaction Caused by Cut-to-Length Forest Operations and Possible Short-Term Natural Rehabilitation of Soil Density. Soil Science Society of America Journal, 75, 2314-2329. https://doi.org/10.2136/sssaj2011.0109
 Manuwa, S. I. (2012). Evaluation of Shear Strength and Cone Penetration Resistance Behavior of Tropical Silt Loam Soil under Uni-Axial Compression. Open Journal of Soil Science, 2, 95-99. https://doi.org/10.4236/ojss.2012.22014
 Moroni, M. T., Carter, P. Q., Ryan, & D. A. J. (2009). Harvesting and Slash Piling Affect Soil Respiration, Soil Temperature, and Soil Moisture Regimes in Newfoundland Boreal Forests. Canadian Journal of Soil Science, 89, 345-355.
 Murphy, P. N. C., Castonguay, M., Ogilvie, J., Nasr, M., Hazlett, P., Bhatti, J. S., & Arp, P. A. (2009b). A Geospatial and Temporal Framework for Modeling Gaseous N and Other N Losses from Forest Soils and Basins, with Application to the Turkey Lakes Watershed Project, in Ontario, Canada. Forest Ecology and Management, 258, 2304-2317.
 Murphy, P. N. C., Ogilvie, J., & Arp, P. A. (2009a). Topographic Modelling of Soil Moisture Conditions: A Comparison and Verification of Two Models. European Journal of Soil Science, 60, 94-109. https://doi.org/10.1111/j.1365-2389.2008.01094.x
 Murphy, P. N. C., Ogilvie, J., Meng, F.-R., White, B., Bhatti, J. S., & Arp, P. A. (2011). Modelling and Mapping Topographic Variations in Forest Soils at High Resolution: A Case Study. Ecological Modelling, 222, 2314-2332.
 Page-Dumroese, D., Jurgensen, M., Elliot, W., Rice, T., Nesser, J., Collins, T., & Meurisse, R. (2000). Soil Quality Standards and Guidelines for Forest Sustainability in Northwestern North America. Forest Ecology and Management, 138, 445-462.
 Roberts, S. D., Harrington, C. A., & Terry, T. A. (2005). Harvest Residue and Competing Vegetation Affect Soil Moisture, Soil Temperature, N Availability, and Douglas-Fir Seedling Growth. Forest Ecology and Management, 205, 333-350.
 Rücknagel, J., Götze, P., Hofmann, B., Christen, O., & Marschall, K. (2013). The Influence of Soil Gravel Content on Compaction Behaviour and Pre-Compression Stress. Geoderma, 209-210, 226-232. https://doi.org/10.1016/j.geoderma.2013.05.030
 Smith, C. W., Johnston, M., & Lorentz, S. (1997). The Effect of Soil Compaction and Soil Physical Properties on the Mechanical Resistance of South African Forestry Soils. Geoderma, 78, 93-111. https://doi.org/10.1016/S0016-7061(97)00029-3
 Soane, B. D., & van Ouwerkerk, C. (1994). Chap.1. Soil Compaction Problems in World Agriculture. In: Developments in Agricultural Engineering (pp. 1-21). Amsterdam: Elsevier. https://doi.org/10.1016/B978-0-444-88286-8.50009-X
 Solgi, A., Naghdi, R., Labelle, E. R., & Zenner, E. K. (2018). The Effects of Using Soil Protective Mats of Varying Compositions and Amounts on the Intensity of Soil Disturbances Caused by Machine Traffic. International Journal of Forest Engineering, 29, 199-207. https://doi.org/10.1080/14942119.2018.1527174
 Sørensen, R., Zinko, U., & Seibert, J. (2006). On the Calculation of the Topographic Wetness Index: Evaluation of Different Methods Based on Field Observations. Hydrology and Earth System Sciences, 10, 101-112. https://doi.org/10.5194/hess-10-101-2006
 Southee, F. M., Treitz, P. M., & Scott, N. A. (2012). Application of Lidar Terrain Surfaces for Soil Moisture Modeling. Photogrammetric Engineering and Remote Sensing, 78, 1241-1251. https://doi.org/10.14358/PERS.78.11.1241
 Vaz, C. M. P., Bassoi, L. H., & Hopmans, J. W. (2001). Contribution of Water Content and Bulk Density to Field Soil Penetration Resistance as Measured by a Combined Cone Penetrometer—TDR Probe. Soil & Tillage Research, 60, 35-42.
 Vega-Nieva, D. J. D., Murphy, P. N. C., Castonguay, M., Ogilvie, J., & Arp, P. A. (2009). A Modular Terrain Model for Daily Variations in Machine-Specific Forest Soil Trafficability. Canadian Journal of Soil Science, 89, 93-109. https://doi.org/10.4141/CJSS06033
 Vogelmann, E. S., Reichert, J. M., Prevedello, J., & Awe, G. O. (2013). Hydro-Physical Processes and Soil Properties Correlated with Origin of Soil Hydrophobicity. Ciência Rural, 43, 1582-1589. https://doi.org/10.1590/S0103-84782013005000107
 Yin, X., & Arp, P. A. (1994). Fog Contributions to the Water Budget of Forested Watersheds in thé Canadian Maritime Provinces: A Generalized Algorithm for Low Elevations. Atmosphere-Ocean, 32, 553-566. https://doi.org/10.1080/07055900.1994.9649512