Analyzing and Projecting Soil Moisture and Cone Penetrability Variations in Forest Soils

This article details how forest soil moisture content (MC) and subsequent re-sistances to cone penetration (referred below as Cone Index, CI) vary by daily weather, season, topography, site and soil properties across eleven harvest blocks in northwestern New Brunswick. The MC- and CI-affecting soil variables refer to density, texture, organic matter content, coarse fragment content, and topographic position (i.e., elevation, and the seasonally affected cartographic depth-to-water (DTW) pattern). The harvest blocks were tran-sect-sampled inside and outside their wood-forwarding tracks at varying times throughout the year. In detail, 61% of the pore-filled moisture content (MC PS ) determinations inside and outside the tracks could be related to topographic position, coarse fragments, bulk density, and forest cover type specifications. In addition, 40% of the CI variations could be related to soil depth, MC PS , and block-specific cover type. Actual versus model-projected uncer-tainties amounted to ΔMC PS ≤ ± 15% and ΔCI ≤ ± 0.5 MPa, 8 times out of 10. Block-centered MC and CI projections were obtained through: 1) daily hydrological modelling using daily precipitation and air temperature weath-er-station records nearest each block, and 2) digitally mapped variations in soil properties, elevation, DTW and forest cover type, done at 10 m resolution.

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

Lowlands (LL)
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.

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

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).
Soil density samples were also collected from each plot by tapping and extracting a metal cup of known volume (85 cm 3 ) 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 D b also containing roots and coarse fragments was determined as per Equation (2).
Weight of dried soil g Weight of Coarse Fragments g D Volume cm The fine-earth bulk density of the soil between rocks and roots was estimated using the oven-dry weight fractions of sand (Sand W ), organic matter (OM W ) and mineral soil depth (cm) as predictor variables as per Equation (3) (Balland et al., 2008).  Since CI tends to be related to pore-filled soil moisture content (MC PS , see Vega-Nieva et al., 2009), it was important to infer soil pore space (PS) and MC PS from the soil particle and bulk densities as follows (Balland et al., 2008): with organic and mineral particle densities set to D OM = 1.3 and D Min = 2.6 g•cm −3 , respectively.

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 (  Figure A1). 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 (FC W ), given by where FC W , SP W , and Sand W 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 (     DTW, elevation, forestcover, soil CI f properties, block, track, brushmat

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 MC V , MC PS and CI as dependent variables. The resulting regression models for Equations (9) , 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.

General Observations and Trends
The plot-based determinations for sand, silt, and clay, OM, CF, MC PS , 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 Table A2 (Appendix). These plots reveal that the MC PS , CI and D b 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 MC PS , D b , 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 MC PS , D b , and CI surveying methodologies, e.g., using D b 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 D b , CI, and MC PS for the top mineral soil 15 cm were somewhat lower and less variable   under matted than non-matted tracks ( Figure 8). Systematic increases for D b 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.  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 MC PS and CI versus CF, sand, OM and D b are presented in Figure 9. As shown, increasing CF generally lowered MC PS 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 MC PS trend with increasing D b relates to decreasing pore space.
Similarly, decreasing MC PS 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 MC PS due to increasing hydrophobicity as OM-containing soils dry out (Vogelmann et al., 2013). In contrast, the increasing MC PS 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 D b (Equation (3); Carter et al., 2007).
Analyzing the trends among the MC PS 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, D b , elevation and tracks, but negatively associated with increasing soil organic matter content.
2) Factor 2 is positively associated with higher elevations where tolerant     3) Factor 3 is positively associated with increasing MC PS 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.

Temporal Soil Moisture Derivations
The ForHyM-generated results for daily soil moisture, stream discharge, snowpack depth and depth of soil frost are illustrated in Figure 11, 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 MC PS0 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 ( Figure A1). These calibrations involved adjusting the For-HyM parameters for snowpack depth and lateral versus vertical layer-by-layer soil permeability, as detailed in Table 5.
Fitting the actual plot-specific MC PS determinations for the top 15 cm of mineral soil versus the corresponding log 10 DTW values with FIA = 1 ha across all the plots by way of Equation (9)     The corresponding best-fitted MC V results, also listed in Table 6, led to:  As to be noted, there is a high to low D b regression coefficient and significance change from MC PS [Equation (14)] to MC V [Equation (15)] due to the greater dependence of MC PS on D b via Equation (5) and Equation (6) while the significance levels for the predictor variables remain about the same. Also note that the best-fitted R 2 and RMSE values drop from MC PS to M V due to the narrower MC V 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 MC PS0 (and therefore MC PS,DWT ) increase, as to be expected.
2) MC decreases towards higher and generally steeper elevations.
3) MC is higher inside than outside soil tracks. The corresponding actual versus best-fitted scatterplots associated with MC PS 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 MC PS 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 MC PS generally higher than inside track.
Grouping the off-and in-track projections and corresponding field determinations into 10% MC PS classes produced the MC PS confusion matrix and associated MC PS cumulative conformance plots in Figure 14. In summary, the ±5% MC PS conformance level so found amounts to 34% across the 20 < MC PS < 100% range, but increases to 81% and 96% as the conformance uncertainty increases to ΔMC PS ≤ ±15 and ±25%, respectively. The possibility of randomly achieving M.-F. Jones, P. Arp Figure 12. Scatterplots of actual versus best-fitted MC PS regression results. Using Equation (12) and (13) demonstrates the extent of block-specific DTW projections on MC PS . 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 MC PS can also be generated from digitally mapped variables (DTW, OM DSM ) 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.

CI Derivation
The best-fitted CI results using plot-determined and plot-projected data are listed in Table 7     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 MC PS 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  (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 log 10 CI values relate to the corresponding log 10 CI projections.   (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 MC PS kappa values, this suggests that the CI projections are slightly more random than the MC PS projections.

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 (MC PS0 ) 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, D b , OM, CF, vegetation type and % canopy closure.
The MC and CI results so produce are projected in Figure 18 for hardwood

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

Conclusion
The above soil MC and CI assessment is limited to forested areas within north-