Relating Cone Penetration and Rutting Resistance to Variations in Forest Soil Properties and Daily Moisture Fluctuations ()
1. Introduction
The soil cone index (CI), a measure of a soil’s resistance to penetration (MPa), is a commonly used soil mechanical property to determine soil strength [1] [2] . This strength generally increases with increasing clay, coarse fragment (CF), and soil density (Db), or reduced pore space (PS), but decreases with increasing soil moisture (MC) and organic matter content (OM, %) [3] [4] [5] [6] . Hence, non-cohesive soils such as sands and sandy loams are more easily penetrated than clay soils [3] [7] [8] , wet soils have low penetration resistances and the resistance to penetration is low for organically enriched soils but high for stony and frozen soils [9] [10] [11] .
In practice, off-road traffic may increase soil compaction and CI, which negatively affects the growth of crops by way of reduced root development [8] [12] [13] [14] [15] . In urban developments, increased CI due to soil compaction decreases soil infiltration of water and tree root growth [16] [17] . However, sufficient CI-index soil strength is needed to allow on- and off-road traffic in agriculture and forestry operations [18] [19] , while off-road recreational traffic needs to be controlled to avoid soil rutting. In this, the resistance of soils to rutting is directly proportional to the ratio between tire footprint pressure and CI [20] [21] [22] . The former increases with increasing vehicle weight and load and decreasing tire footprint, which―in turn―decreases with increasing tire width, wheel diameter, and decreasing tire pressure. In the field, rut depths further increase from single to multiple passes, and with slope-in- duced tire spinning [23] .
Efforts to minimize soil rutting require reliable forecasting of off-road soil trafficability. Doing this, however, is challenging because soil and machine-use conditions may vary daily from location to location. By location, low CI conditions do not last as long for sandy soils than for loams and clays. In addition, soil trafficability varies with the extent of soil freezing and thawing, especially when traffic turns thawing soils into mud [24] .
The objective of this article is determining how manually derived soil CI determinations change in response to weekly spring-to-winter changes in soil moisture and temperature for three contrasting soil conditions. The data so generated allowed for: 1) quantifying the relationship between CI and soil MC; 2) emulating and interpreting the changes in soil moisture, CI, and rutting depth; 3) daily year-round modelling of soil trafficability by soil texture and soil depth. While machine-based cone penetration testing (CPT [25] ) would be more accurate and precise, manual CI determinations have the greater portability and affordability advantage for assessing how soil trafficability conditions vary from location to location across landscapes and seasons.
2. Materials and Methods
2.1. Location Description
Three forest sites in Fredericton, New Brunswick, were chosen for this study (Figure 1, Table 1):
1) A mixed-wood stand on sandy clay loam in a wooded section on the University of New Brunswick campus (UNB);
2) A hemlock (Tsuga canadensis) stand on a rich loam in Odell Park (OP);
3) A silver maple site (Acer saccharinum) on an alluvial sandy loam next to a fresh-water marsh within the floodplain of the Nashwaaksis stream (SM).
The two non-alluvial soils developed on grey sandstone ablation / basal till. Elevation for the three sites ranges from 6 to 70 m [26] . The topography varies from undulating to hilly. The upland forest vegetation is representative of the Acadian forest species, i.e., sugar maple (Acer saccharum), red maple (Acre rubrum), white birch (Betula papyrifera), balsam fir (Abies balsamea), black spruce (Picea mariana), and hemlock (Tsuga canadensis). The 1950-2017 Fredericton weather record has a mean annual temperature of 6.6˚C, with monthly means of −1.8 and 14.9˚C for January and July, respectively. Mean annual precipitation amounts to 1100 mm, including 250 mm of snow [27] .
2.2. Field Experiment
Soil layers were described and samples were taken from two freshly dug soil pits at each of the three locations. Five soil volumetric moisture content (MCy) and CI readings were taken manually each week from May 29, 2015 to November 2, 2015 within two circular plots (1.5 m radius) near the soil pits at each location. This was done using a Delta T HH2 moisture meter and a Humboldt digital cone penetrometer (cone are at base = 1.5 cm2; cone angle 60˚). The MCv readings were taken at 7.5-cm mineral soil depth. The CI readings were obtained at 15, 30, 45, and 60 cm depths, but were not recorded where obstructed by logs, coarse roots, and surface-accumulated rocks.
![]()
Figure 1. Overview depicting of the three Fredericton (New Brunswick) locations (left), and site-specific plot locations for SM (top), OP (middle), UNB (bottom) (Imagery Source: Esri, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community).
![]()
Table 1. Location descriptions used for initializing ForHyM.
The soil samples were placed into labeled freezer bags for storage. Prior to analysis, the samples were dried in a forced-air oven 75˚C for 24 hours, crushed with a mortar and pestle, and passed through a 2-mm sieve to remove and to determine the CF. The fine-earth fraction was used to determine its sand, silt, and clay content using the hydrometer method [28] . The soil carbon content (C) of this fraction was determined using a LECO CNS-2000 analyzer. Soil OM content was estimated by weight by setting OMg% = 1.72 × C%. The pore-space filled moisture content (MCps) was inferred by assuming that soil gravimetric moisture content (MCg), soil bulk density (Db) and the PS percentage would be affected by depth and OM content as follows [3] :
(1)
(2)
(3)
where Dp is particle density (2.65 g/cm3), and PS is the pore space fraction of the fine earth.
2.3. Hydrological Modelling
The forest hydrology model (ForHyM) [29] [30] [31] was used to emulate the changes in daily soil moisture, soil temperature and snowpack conditions for each of the three locations from 2006 to 2017. Doing this involved compiling the daily Fredericton weather records for air temperature, precipitation (rain, snow), stream discharge, and open-ground snow depth [27] [32] . Also specified were elevation, slope, aspect, and extent of forest cover (Table 1). The model-internal water and heat flow parameters pertaining to soil permeability, thermal conductivity, and heat capacity were plot-adjusted by texture, OM and CF content (Table 2), and by comparing actual with modeled soil moisture content. This was done through manually resetting the default values for: 1) the air-to-snow- pack heat-transfer coefficient; 2) the initial snowpack density of freshly fallen snow to reflect the open-ground conditions at the weather station [33] ; and 3) the lateral soil permeability to account for lateral flow tortuosity [34] [35] . These adjustments ensured that the model output conformed to actual snowpack depth and stream discharge records.
2.4. Data Analysis and Model Projections (MCv, CI, Rut Depth)
The data and ForHyM estimates for MCv, CI, texture, CF, OM, Db, and PS were entered into a spreadsheet by location, date, and soil depth. This compilation served 1) to generate basic statistical summaries, 2) to analyze the measured and modelled time-series plots for MCv and CI, and 3) to determine the best-fitted linear and multiple regression models with CI as dependent variable, and with MCv (measured, modelled), soil texture, OM, CF, PS, and soil depth as independent variables. A linear regression model served to relate measured CI at 15, 30, 45 and 60 cm soil depth to measured and ForHyM-modelled MCv. A multiple regression model served to relate CI to MCv, PS, and CF as follows:
(4)
where MCps is the water-filled portion of the PS, in percent. The best-fitted model so generated was incorporated into the ForHyM model to determine how MCy, CI and rutting depths pertaining to all-terrain vehicle (ATV) traffic would vary over time at each of the three locations. The equations adopted for rut modelling were as follows [4] [36] :
Potential rut depths for n passes:
(5)
with NCI (the nominal cone index) given by:
(6)
where b is tire width (m), d is tire diameter (m), h is section height (m), δ is tire deflection (m) given by 0.008 + 0.001 (0.365 +170 / p ), p is tire inflation pressure (kPa), W is vehicle weight + load (kN) per wheel, and n is number of vehicle passes along the same track. Potential rutting depths for all-terrain recreational vehicle (ATV) traffic were determined using the following machine specifications: number of wheels = 4; W per wheel = 3.1 kN; b = 0.254 m; d = 0.62 m; h = 0.3 m; p = 34.4 kPa; n = 10 passes.
![]()
Table 2. ForHyM initialization requirements by soil layer per plot and location.
To visually represent the temporal changes in MC topographically and over seasons, MCps was spatially related to the depth-to-water index (DTW). This as index was generated from a 1-m resolution bare-earth digital elevation model (DEM) for the Fredericton area [37] . This index determines the elevation rise along the least slope path from each cell across the landscape to its nearest open-water cell corresponding to streams, lakes, rivers and open shores [38] [39] . Changing the upslope flow-accumulation area by channel flow initiation (FI), i.e., changing the amount of upstream area needed to initiate streamflow, allows for indexing DTW by season. For example, FI = 4 ha generally represents permanent stream flow at the end of summer, FI = 0.25 ha represents the extent of ephemeral stream flow during and after snowmelt, and FI = 1 ha represents channel flow during the transitional periods from fall to winter. The resulting DTW rasters with FI = 4, 1 and 0.25 ha were used to determine how the soil moisture conditions and rutting depths would vary across the terrain associated for the three sampling locations by season. This was done by applying Equation (7) and Equation (8) [3] , i.e. [40] :
(7)
(8)
with p = 2, p as soil-specific parameter ranging from 0.2 to 2 and), and DTWridge (in m) RDn,ridge and RDn,DTW=0 (in mm) determined for driest and wettest parts of each.
3. Results
3.1. Soil Moisture and CI Measurements
Each of the three locations showed distinct variations in soil properties, strength, and moisture readings over the course of 23 weeks. Given the plot-by-plot soil property differences―and tracking the changes in soil moisture over time―re- vealed that the OP plots drained quickly. In contrast, the UNB plots varied the most from wet to dry and back again to wet from spring to fall (Figure 2). In direct correspondence, resistance to cone penetration varied the least for the two SM plots, and the most for the UNB plots. These differences arose from the compacted and poorly drained sandy loam for the UNB plots, the well-drained loamy sand with low CF content for the OP plots, and seasonally recurring flooding of the SM plots (Table 1 and Table 2). The high springtime levels for MCv within the top 15-cm soil at the UNB and SM locations are due to high Ah-layer OM content, which―according to Equation (1)―lowers Db and enhances the soil-filled PS between the coarse fragments.
Plotting the CI measurements at 15 cm depth to the MCv measurements revealed that the log-transformed CI and MCv values are linearly related to one another as shown in Figure 3 (left), as follows:
(9)
R2 = 0.60, RSME = 0.13, MAE = 0.10, with the UNB location coded 1 and 0 otherwise. Similarly similar strong correlations between MC and CI have been reported elsewhere [5] [41] [42] . With respect to increasing soil depth―and as
![]()
Figure 2. Left: Measured MCv for the top 15 cm of soil. Right: Measured CI at 15, 30, 45, and 60 cm depth for Plots 1 and 2 at the OP, UNB, and SM locations.
![]()
Figure 3. Scatterplots of measured log10CI versus measured log10MCv (left), and weekly of CI / CImax averages for the OP, UNB, and SM locations (right).
shown in Figure 3 (right)―CI increases, by plotting the ratio of the weekly averages of CI over CImax per plot by location. A similar trend has also been reported elsewhere [43] [44] [45] .
3.2. Soil Moisture and CI through Hydrological Modelling
The modelling of the year-round soil moisture conditions required plot-specific ForHyM initializations and calibrations. These included the Fredericton-specific calibrations for snowpack depth and stream discharge required using daily Fredericton Airport weather records for rain, snow and air temperature, and adjusting the ForHyM-default settings for lateral and downward water flow, as listed in Table 3. The plot initializations in Table 1 and Table 2 refer to entering the plot- and/or layer-specific values for slope, aspect, vegetation type and cover, forest floor depth, percentages for sand, silt, clay, CF, OM, and layer depth.
Shown in Figure 4 are the resulting time-series plots for daily air temperature and precipitation (input), snowpack depth, stream discharge, top 15-cm soil MCv (actual and modelled), and frost depth (modelled). The resulting scatter plots in Figure 5 for actual and best-fitted ForHyM snowpack depth and top 15-cm MCv (top 15 cm) demonstrate a reasonable good fit, with R2 = 0.81 for snowpack depth, 0.62 for stream discharge, and 0.76 for MCv (Table 4).
For the purpose of predicting how CI would vary across time by soil texture, Db and CF content (Table 2), it was necessary to use the ForHyM-generated depth- and time-dependent MCv output for the 0 - 15, 15 - 30, 30 - 45 and 45 - 60 cm soil layers as predictor variable. Doing this involved estimating how much of the infiltrating and percolating water would be retained at any time within the fine-earth fraction between the coarse fragments of each layer. For example, the space available for water retention would decrease with increasing CF content. Consequently, there would be less PS to fill between the coarse fragments during wet weather conditions, and there would also be less water available for root uptake during warm summer weather [46] . This being so, The ForHyM-generated projections in Figure 6 by location and soil layer show greater MCv and MCps variations for the stony UNB location, followed by the less stony SM and the more sandy OP locations. In combination, the ForHyM projections in Figure 6 capture the plot-by-plot MCps variations such that OPMC > SMMC > UNBMC.
Figure 7 and the correlation coefficients in Table 5 show how CI varies with varying soil texture (Sand), CF, OM, PS and MCps. In general, CI decreases with
![]()
Table 3. ForHyM calibrations for the Fredericton area: default multipliers.
![]()
Figure 4. ForHyM time-series plots for daily air temperature and precipitation (ForHyM input), actual and modelled output for stream discharge and snowpack depth, and location-specific modelled frost depth (modelled).
increasing PS and sand content due decreasing particle-to-particle contacts. Increased OM content decreases CI by way of soil aggregation, i.e. by further loosening the point of contact among the aggregated soil particles. The CF-induced increase on CI refers to the increasing strength needed to displace the coarser particles away from cone penetration path [47] . Together, Sand, OM, CF and PS affect the daily variations in CI and soil moisture retention through their combined effect on soil pore space, texture, structure and drainage [33] [48] [49] .
Subjecting the correlation matrix in Table 6 to factor analysis revealed that the CI variations can be grouped into three CI-determining factors. Factor 1 is the Location Factor, which relates a component of the CI variations to the location- and layer-specific CF and PS determinations. Factor 2 is the Soil Moisture Factor, which relates some of the CI variations to MCps. Factor 3 is strongly related to Sand, but―in this formulation―has no salient effect on CI.
![]()
(a)![]()
![]()
(b) (c)
Figure 5. Actual versus ForHyM best-fitted scatter plots for MCv (top 15 cm) (left), monthly stream discharge (middle), and daily snowpack depth (right).
![]()
Table 4. Best-fitted regression model for measured (actual) versus modelled top 15-cm soil MCv by location (UNB, SM, OP) and overall.
Using PS, MCps, and CF as independent variables produced the following best-fitted multiple regression result for all soil layers and locations combined:
![]()
(10)
R2 = 0.54, RMSE = 0.36, MAE = 0.29. This result is illustrated in Figure 8 by way of the 3D plots, which reveal moderate CI increase with decreasing MCps, and a rapid CI increase with increasing CF. In reality, CI and soil strength should decrease again as MCps drop towards zero as the soil becomes more brittle due to reduced particle-to-particle hydrogen-bonding at low MC [50] .
![]()
Figure 6. ForHyM-generated MCv and MCps projection for the 0 - 15, 25 - 30, 30 - 45 and 45 - 60 cm soil layers by plot and location.
![]()
Figure 7. Plotting plot-by-plot measured CI vs. OM, Sand, PS, MCps, and CF, showing PS, MCps, and CF as stronger CI predictor variables than OM and Sand.
![]()
Table 5. Correlation matrix for plot- and layer-determined CI, OM, Sand, CF and ForHyM-estimated Db, MCv, MCps.
![]()
Table 6. Factor analysis of Table 5.
![]()
![]()
(a) (b)
Figure 8. Modelled CI (Equation (10)) in relation (a) to MCps and PS at CF = 20%; and (b) to CF and MCps at SP = 20%.
While Sand and OM are important water retention and porosity predictor variables [51] [52] , including them as part of the multiple regression process did not significantly improve the best-fitted results, likely due to the significant correlations between OM and PS and between Sand and CF in Table 5. However, adding sampling location to the predictor variables (each location coded 1 where applicable, else 0) improved the best-fitted result as follows:
![]()
(11)
R2 = 0.60, RMSE = 0.33, MAE = 0.27. This means that the CI values at the SM plots are, on average, slightly lower than at the other locations. This difference may be related to unaccounted differences pertaining to, e.g., CF size (generally smaller at SM than at the other two locations), and differences in rooting pattern.
Repeating this analysis by location and by soil depth produced the best-fitted results listed in Table 7. From this, it can be noted that R2 remained about the same by location, varying from 0.41 (SM) to 0.66 (UNB), but decreased with increasing soil depth from 0.68 at the top to 0.10 at 60 cm soil depth. This decrease would mostly be due to the location-by-location Db, MC and CF differences. This is because 1) the ForHyM-generated MC estimates already take the effect of CF on MCps into account, and 2) the CI readings become increasingly erratic when pushed through soils with increasing CF content.
The dependency of CI data on soil PS, MC and CF content was further evaluated through multiple regression analysis based on literature-generated CI formulations (Table 8). The result of so doing indicated that: 1) Equation (10) provides the best data representation overall, 2) the linear formulations for CI are somewhat weaker than the logarithmic formulations. Also, 3) soil porosity (or density) and MC are the more persistent and significant CI predictor variables than either Sand or CF alone.
3.3. Predicting Potential ATV-Caused Soil Rutting Depth
ForHyM was used to transform the MCps and CI projections over time into likely ATV-generated rut depths over time from April 2013 to April 2017, using the average top 15-cm PS and CF values and Equation (7), Equation (8), and Equation (10) the two plots at the three sampling locations. The results are represented by the time-series plot in Figure 9. As to be expected, deepest ruts would be incurred during spring and fall, with minor blips during summer. Ruts could also be incurred during winter when some of the frozen soils would thaw due to interim warm weather and upward geothermal heat flow underneath the heat?insulating snow accumulations [53] . While trafficability advisories exist from fall to spring due to wet soil conditions, such advisories apply regionally, and therefore fall short in terms of local “when” and “where” decisions.
![]()
Table 7. Linear regression results for measured vs. modelled CI by depth and location.
![]()
Table 8. Review of functional relationship between CI and soil properties.
1: This study; 2: [4] ; 3: [45] ; 4: [8] ; 5: [54] ; 6: [55] .
The extent to which soil rutting would be seasonally affected across the general neighbourhood of each of the three location was ascertained through digitally generating the elevation-derived cartographic depth-to-water index (DTW) associated with the 4, 1 and 0.25 ha upslope areas for streamflow initiation [40] (Figure 10). Using these patterns in combination with Equations (7) and Equations (8) produced the spatial MCps and potential ATV-related rut-depth maps in Figure 11, intended to be representative of the off-road soil trafficability conditions during spring, end of summer and the fall to winter transition. As shown, the UNB location has the potential to be the most trafficable among the three locations in summer, but would be worst during spring and fall. In contrast, the OP location would have the least traffic impact across the area and seasons based on texture-facilitated soil drainage. However, moderate soil rutting could occur within the 4-ha DTW < 1 m zone at OP. Overall, the soil rutting conditions follow these sequences: dry weather: UNB < OP < SM; wet weather: OP < SM < UNB (Figure 9).
4. Discussion
This article describes ways and means by which the resistance of soils to cone penetration can be analyzed and modeled at the daily level year-round, over many years, and for the varying soil conditions by select locations. The results so obtained are―apart from study-specific biases―generally consistent with what has been reported in the literature. These biases would inter alia refer to differences in CI methodology by, e.g., cone dimensions, speed of cone penetration, and field versus laboratory testing [4] .
While the plot-by-plot determinations of this study are limited to three contrasting forest locations, they are at least representative of how soil moisture, CI, and rutting depth vary by soil properties, season and topographic position, as demonstrated through daily and spatial modelling. The extent to which this
![]()
Figure 9. ForHyM-generated unfrozen MCps, CI, and rut depths from April 2013 to April 2017 for the topsoil (top 15 cm of soil) for plot 1 (top) and plot 2 (bottom) at UNB, OP, and SM.
approach can be generalized requires additional research. For example, the spatial and DTW-dependent soil trafficability formulation for CI and rut depth should be tested across a wider range of glaciated and non-glaciated landforms. Doing so would involve extending the above regression analyses across a wider range of independently varying soil types and properties. For example, where soils are cemented because of pedogenic Fe and Ca accumulations, the approach
![]()
Figure 10. Cartographic depth-to-water index (DTW ≤ 1 m), overlain on the hill-shaded LiDAR-derived bare-earth digital elevation model for the UNB, OP, and SM locations for end-of-summer (top), spring-to-summer as well as fall-to-winter (middle), and early- spring, as emulated using upslope stream-flow initiation areas amounting to 4, 1 and 0.25 ha, respectively.
would need a cementation predictor variable. In some cases, the mix of the best-fitting regression variable and regression coefficients may also differ, as demonstrated above in Table 8.
Key to applying the approach across time and landscapes is the ability to estimate how soil trafficability changes in direct response to the spatially and temporally varying topo-pedo-hydrological conditions, meter-by-meter. Traditional soil survey maps can be helpful in this regard but only if the individual map units and borders conform to actual soil drainage contours. To this extent, further progress can be made by:
1) refining and adjusting each unit to its landform- and DEM-defining drainage position;
2) exploring how the trafficability affecting soil properties (MC, texture, CF, OM, Db, depth) vary across the landscape of interest from the highest to the lowest elevation points;
3) determining the point of streamflow initiation inside each flow channel either through field observations or through DEM-based flow-initiation algorithms.
![]()
Figure 11. Soil moisture content per pore space [MCps (%)] and all-terrain rut depth after 10 passes along same track (RD10, mm), generated from the season?representative DTW patterns in Figure 10 using Equation (7) and Equation (8). Top: end-of summer. Middle: spring-to-summer and fall-to-winter transitions. Bottom: after snowmelt.
Together, these refinements would add further precision to the soil moisture and rut depth maps in Figure 11. For example, there would be a noticeable difference between DTW, MCps and ATV rut depth projections within and outside the floodplain associated with the SM location.
Some progress towards these refinements has already been made in terms of checking existing trail conditions in terms of ATV-induced rutting extent, and by correlating this extent to the ridge-to-valley of the cartographic depth-to- water index (DTW [40] ). The multi-pass implications on wood-forwarding rutting depth have been reported by [56] , and were further evaluated by [4] by way of Equation (7) and Equation (8). However, much more work needs to be done by not only addressing the DTW-emulated variations in soil wetness but also by addressing the changes in Db, texture, CF and OM content as these would vary from ridge tops to valleys in a systematic manner. For example, upslope soils would generally be thinner and coarser with less OM than downslope soils. The reverse would occur in severely eroded medium-textured soils, with the more cohesive soil remains upslope and the more easily eroding sand and silt fractions accumulating downslope.
Since the above analysis is restricted to bare ground conditions and mineral soil layers, rut-reducing surface accumulations of snow, ice, forest litter, peat, and roots are not addressed. Bare-ground conditions, however, exist across forested landscapes along non-paved roads, after ground-exposing operations such as root extractions, mounding and plowing, and underneath forest cover where litter accumulations are low or absent due to fast litter decomposition rates. The latter condition is more prevalent under hardwood and pine forests than under fir and spruce forests. Repeated recreational traffic in such areas under moist to wet weather conditions would induce significant rut-induced damage through trail braiding, soil erosion, gulley formation, and stream and lake sedimentation [45] .
Also not addressed are the effects of snow and ice build-up on top of soils during winter, which would increase the resistance to soil penetration, compaction, and rutting through increased load-bearing capacities. Since not all the water is frozen in sub-zero clay- and OM-enriched soils, there could be problems associated winter-based soil rutting followed by instantaneous flash freezing. In summary, the above soil rutting assessment is only applicable for bare ground conditions. Soils covered by forest litter, slash, snow, and ice would obviously reduce rutting.
5. Concluding Remarks
The above soil rutting assessment via manual testing of the temporal changes in the soil resistance to penetration is limited to the immediate area at and around the three sampling locations of this study. More research is needed to extend and test this research regarding general applicability. As shown, the approach taken would allow this by way of hydrological and digital elevation modelling, and further procurement of CI-relevant soil information.
Acknowledgements
This research was supported by a Collaborative Research Developement Project
sponsored by J.D. Irving, Limited, and the Natural Science and Engineering Council of Canada. Special thanks go to Doug Hiltz, Friedrich Wüthrich, and Brittany Hartery for weekly field sampling.