Relating Cone Penetration and Rutting Resistance to Variations in Forest Soil Properties and Daily Moisture Fluctuations

Soil resistance to penetration and rutting depends on variations in soil texture, density and weather-affected changes in moisture content. It is therefore difficult to know when and where off-road traffic could lead to rutting-induced soil disturbances. To establish some of the empirical means needed to enable the “when” and “where” determinations, an effort was made to model the soil resistance to penetration over time for three contrasting forest locations in Fredericton, New Brunswick: a loam and a clay loam on ablation/ basal till, and a sandy loam on alluvium. Measurements were taken manually with a soil moisture probe and a cone penetrometer from spring to fall at weekly intervals. Soil moisture was measured at 7.5 cm soil depth, and modelled at 15, 30, 45 and 60 cm depth using the Forest Hydrology Model (ForHyM). Cone penetration in the form of the cone index (CI) was determined at the same depths. These determinations were not only correlated with measured soil moisture but were also affected by soil density (or pore space), texture, and coarse fragment and organic matter content (R2 = 0.54; all locations and soil depths). The resulting regression-derived CI model was used to emulate how CI would generally change at each of the three locations based on daily weather records for rain, snow, and air temperature. This was done through location-initialized and calibrated hydrological and geospatial modelling. For practical interpretation purposes, the resulting CI projections were transformed into rut-depth estimates regarding multi-pass off-road all-terrain vehicle traffic.


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 (D b ), 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-induced 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.

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

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 (MC y ) 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 cm 2 ; cone angle 60˚).The MC v 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.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 OM g % = 1.72 × C%.The pore-space filled moisture content (MC ps ) was inferred by assuming that soil gravimetric moisture content (MC g ), soil bulk density (D b ) and the PS percentage would be affected by depth and OM content as follows [3]: where D p is particle density (2.65 g/cm 3 ), and PS is the pore space fraction of the fine earth.

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

Data Analysis and Model Projections (MCv, CI, Rut Depth)
The where MC ps 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 MC y , 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: with NCI (the nominal cone index) given by: 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  To visually represent the temporal changes in MC topographically and over seasons, MC ps 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]: with p = 2, p as soil-specific parameter ranging from 0.2 to 2 and), and DTW ridge (in m) RD n,ridge and RD n,DTW=0 (in mm) determined for driest and wettest parts of each.

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-revealed 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 MC v 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 D b and enhances the soil-filled PS between the coarse fragments.
Plotting the CI measurements at 15 cm depth to the MC v measurements revealed that the log-transformed CI and MC v values are linearly related to one another as shown in Figure 3 (left), as follows: R 2 = 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

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 MC v (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 MC v (top 15 cm) demonstrate a reasonable good fit, with R 2 = 0.81 for snowpack depth, 0.62 for stream discharge, and 0.76 for MC v (Table 4).
For the purpose of predicting how CI would vary across time by soil texture, D b and CF content (Table 2), it was necessary to use the ForHyM-generated depth-and time-dependent MC v 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 MC v and MC ps 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 MC ps variations such that OP MC > SM MC > UNB MC .
Figure 7 and the correlation coefficients in Table 5 show how CI varies with varying soil texture (Sand), CF, OM, PS and MC ps .In general, CI decreases with Table 3. ForHyM calibrations for the Fredericton area: default multipliers.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 MC ps .Factor 3 is strongly related to Sand, but-in this formulation-has no salient effect on CI.Using PS, MC ps , and CF as independent variables produced the following best-fitted multiple regression result for all soil layers and locations combined: R 2 = 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 MC ps , and a rapid CI increase with increasing CF.In reality, CI and soil strength should decrease again as MC ps drop towards zero as the soil becomes more brittle due to reduced particle-to-particle hydrogen-bonding at low MC [50].R 2 = 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 R 2 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 D b , MC and CF differences.This is because 1) the ForHyM-generated MC estimates already take the effect of CF on MC ps 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.

Predicting Potential ATV-Caused Soil Rutting Depth
ForHyM was used to transform the MC ps 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.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 MC ps 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).

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, 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;  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, MC ps 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-towater 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 D b , 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.

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.
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.4kPa; n = 10 passes.

Figure 2 .
Figure 2. Left: Measured MC v 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 .
Figure 3. Scatterplots of measured log 10 CI versus measured log 10 MC v (left), and weekly of CI / CI max averages for the OP, UNB, and SM locations (right).

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

Figure 7 .
Figure 7. Plotting plot-by-plot measured CI vs. OM, Sand, PS, MC ps , and CF, showing PS, MC ps , and CF as stronger CI predictor variables than OM and Sand.

Figure 9 .Figure 10 .
Figure 9. ForHyM-generated unfrozen MC ps , 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.

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 .
Figure 11.Soil moisture content per pore space [MC ps (%)] and all-terrain rut depth after 10 passes along same track (RD 10 , 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.

Table 1 .
Location descriptions used for initializing ForHyM.
data and ForHyM estimates for MC v , CI, texture, CF, OM, D b , and PS were v (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 MC v .A multiple regression model served to relate CI to MC v , PS, and CF as follows:

Table 2 .
ForHyM initialization requirements by soil layer per plot and location.

Table 4 .
Best-fitted regression model for measured (actual) versus modelled top 15-cm soil MC v by location (UNB, SM, OP) and overall.

Table 5 .
Correlation matrix for plot-and layer-determined CI, OM, Sand, CF and ForHyM-estimated D b , MC v , MC ps .

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.