Soil compaction due to heavy forest traffic: measurements and simulations using an analytical soil compaction model

Abstract

• Context

Models for predictions of soil compaction following forest traffic represent important decision tools for forest managers in order to choose the best management practices for preserving soil physical quality. In agricultural soil compaction research, analytical models are widely used for this purpose.

• Aims

Our objective was to assess the ability of an analytical model to predict forest soil compaction under forwarder traffic.

• Methods

We used the results from two experimental sites set up in north-eastern France in 2007 and 2008 to compare simulations using the SoilFlex model with observed bulk density following forwarder traffic.

• Results

The best model-based predictions were found when considering the mean initial soil conditions and an increased rebound parameter in the upper soil layers (0–10 cm) in comparison to the deeper layers (10–50 cm). The need to increase the rebound parameter in the soil surface layer to improve model accuracy was attributed to a large soil organic matter content in the uppermost layers of forest soils. For the site where initial soil mechanical parameters were measured as a function of soil bulk density and water content, the model performance was good, with a root mean square error (RMSE) of 0.06. The model performed poorer (RMSE of 0.11), especially for the surface soil layer, for the second site that was wetter at the time of traffic and where soil mechanical properties were not measured but estimated by means of pedo-transfer functions.

• Conclusions

SoilFlex was found to yield satisfactory predictions and could help forest managers estimate the risk of compaction and to select the most appropriate machinery for given soil conditions in order to preserve the soil from physical degradation during traffic in forest ecosystems. However, our results emphasise the need for research on soil mechanical properties of forest soils, in particular on the role of soil organic matter and roots on soil compressive properties.

1 Introduction

Nowadays, an increase in mechanisation of forest operations in France threatens soil sustainability, despite increasing awareness and concern. Indeed, heavy forest traffic has been proven to cause soil physical degradation, resulting in changes in soil chemical and biological properties (Greacen and Sands 1980; Horn et al. 2007). In disturbed forest ecosystems, soil structure regeneration occurs through natural process as soil remediation through tillage and/or fertilisation (as in agriculture) are missing in management practices. The natural recovery rates to undisturbed levels of structure and porosity are generally small, thus posing the question of soil resilience after compaction. Therefore, the prevention of soil physical degradation is particularly needed in forest ecosystems. Soil compaction simulation models (Défossez and Richard 2002) are important decision tools for prevention of excessive soil compaction by heavy traffic.

Forest soils differ from agricultural soils by their organic matter content and quality, porosity, acidity, root and gravel content as well as heterogeneity. Therefore, different mechanical behaviours may be expected. Root networks increase soil resistance to deformation under heavy load (Soane 1990; Cofie et al. 2000). However, assessment of several parameters not easily accessible to forest managers (e.g. rooting pattern, root stiffness and strength) is needed to take into account the soil reinforcement effect of roots in compaction modelling (Cofie et al. 2000). The presence of forest floor litter, the increased soil organic carbon (SOC) content and the modified organic matter quality may decrease the compactibility of forest soils in comparison to agricultural soils (Soane 1990). However, the effect of SOC on soil compressibility has been shown to depend on soil water content (Soane 1990; Pereira et al. 2007). Blanco-Canqui et al. (2005) found smaller soil strength and larger SOC content in forest soils than in pasture and cultivated soils. They also observed that increasing SOC increased water retention and decreased soil strength, both increasing soil susceptibility to compaction. Yet, Mosaddeghi et al. (2000) and Williamson and Neilsen (2000) observed that an increase in SOC content lead to an increase in critical moisture content, i.e. the soil water content at which the maximum deformation is achieved for a given stress. Mosaddeghi et al. (2000) stated that increasing SOC content allowed the soil to better withstand loading especially at large water content. However, Pereira et al. (2007) found the opposite. Schäffer et al. (2008) found that the impact of heavy agricultural machinery on soil structure was the same whatever the SOC content. This last result implies that even if increasing SOC could increase the susceptibility to compaction, through an increase in water content and a decrease in soil strength, it did not increase soil deformation. As for forest compaction literature, studies on the role of SOC on soil mechanical properties are scarce and do not allow to draw any clear conclusion or tendency. Vega-Nieva et al. (2009) observed that the effect of soil organic matter on forest soil resistance to penetration (measured in situ) was only due to its effect on soil bulk density. Smith et al. (1997a) found a significant effect of loss-on-ignition (combining the effects of soil texture and organic matter) on compactibility and compressibility. Clear additional effect of SOC besides the effect of soil texture were only observed by the previous authors for soils with small clay contents (<25 %). Smith et al. (1997b) found for South African forestry soils a differential role of soil water content on the compaction susceptibility as a function of SOC and texture. Evidences of the role of SOC on soil susceptibility to compaction arise also in studies on forest soils, yet no relationships have been calibrated allowing to predict soil mechanical properties as a function of SOC and other soil properties. Besides studies on the relationship between soil strength and SOC do not allow to draw any conclusions concerning the plastic behaviour of soils with different SOC content and quality.

Keller et al. (2007) developed a model named SoilFlex that enables calculations of soil stresses and compaction due to agricultural field traffic. SoilFlex proposes different analytical procedures to calculate the soil–tyre contact area, the normal stress distribution at this interface, the vertical stress propagation through the soil profile and the resulting soil deformation (soil displacement, rut depth and bulk density). It is easy to use, requires a limited amount of parameters and also includes several pedo-transfer functions to estimate soil mechanical parameters from easily collectable soil parameters like bulk density, texture and water content. Models based on analytical procedures have been validated for agricultural soils (Défossez and Richard 2002; Keller et al. 2007). Therefore, SoilFlex simulations could help forest managers to choose the harvesting equipment that is the least likely to cause damage to the soil. However, the model requires validation for forest soils, especially as it does not take into account the effects of gravels, roots, and organic matter on stress–strain relationships.

Our objectives were:

• 1. To simulate soil compaction in forest soils due to forwarder traffic by means of an analytical soil compaction model, the SoilFlex model.

• 2. To analyse the effect of variability in soil water content and bulk density on soil compaction simulations.

• 3. To evaluate the impact of organic carbon on bulk density values after wheeling.

2 Materials and methods

2.1 Study sites

Two experimental sites have been set up in Lorraine (NE part of France) to study the effects of soil compaction. They are located in the ‘Hauts-Bois’ forest-Azerailles (AZ; 48°29′19″ N, 6°41′43″ E), Meurthe et Moselle, and in the ‘Grand Pays’ forest-Clermont en Argonne (CA; 49°06′23″ N, 5°04′18″ E), Meuse. The sites of CA and AZ have an elevation of 270 and 300 m, respectively. The climate of the region is characterised by a 30-year mean annual temperature of 9 (AZ) to 9.5 °C (CA) and a 30-year annual precipitation of 900 (AZ) to 1,000 mm (CA).

Selected soil properties of the study sites are presented in Table 1. The soil of both sites is classified as a neoluvisol (ruptic) (IUSS Working Group WRB, 2006) and is developed on a silt loam layer approximately 50 cm thick laying on a clayey layer (weathering of a Keuper marl for the AZ site and weathering of gaize rock for the CA site). This strong textural discontinuity caused limited temporary water logging; at the lower part of the silt, loam layer hydromorphic features could be observed.

Table 1 Selected soil characteristics (mean and standard error)

2.2 Experimental treatment

Each site was clear-cut and timber (Fagus sylvatica and Quercus petraea mainly) was extracted using a cable yarding system to avoid soil damage. All remaining slash was removed by hand or using an iron horse to avoid soil compaction and/or sinking of the slash into the soil during traffic. Afterwards, an eight-wheel drive forwarder (1996 Valmet 840, serial number 9146, Valmet Logging, Sweden) drove on the soil for an equivalent of two passes (one forward and one rearward pass), in May 2007 at AZ and in March 2008 at CA. The wheel tracks were adjacent to each other in order to create an equally levelled compacted area of 30 × 50 m. In each of the three blocks designed to account for spatial heterogeneity, land strips of the same surface area (30 × 50 m) adjacent to the trafficked plots (T) remained undisturbed by the forwarder and were considered as control plots.

The tyres of the forwarder were 60 cm wide, had a diameter of 133 cm (600/55 × 26.5) and were inflated to a pressure of 360 kPa for both sites. The empty forwarder weighed 11.4 Mg, the four front wheels supporting 6.9 Mg and the four rear wheels supporting 4.5 Mg. In AZ, the wood-loaded forwarder weighed 23.3 Mg, the four front and the four rear wheels supporting 7.56 Mg (i.e. the empty weight on the four front wheels + 5 % of the load) and 15.76 Mg (i.e. the unloaded weight on the four rear wheels + 95 % of the load), respectively. In CA, we only weighed the wood load and deduced the total weight of the loaded forwarder (16.7 Mg), the loaded weight on the four front wheels (7.17 Mg) and rear wheels (9.57 Mg) according to the measurements taken in AZ.

2.3 Measurements of soil compaction

2.3.1 Initial state

For each soil layer, we measured soil gravimetric water content (WC in g 100 g⁻¹) and bulk density (BD in grammes per cubic centimetre), which were used to obtain the soil mechanical properties using the pedotransfer functions developed by Saffih-Hdadi et al. (2009) (see below).

Before setting up the experimental sites (in 2006 for both sites, under previous forest cover), 18 and 13 pits at AZ and CA, respectively, were sampled to determine BD using steel cylinders (10 cm height and 5 cm diameter) at 0–145 cm depth using 10 cm depth increment. These pits were spread over the entire site area that was not disturbed by the following forwarder traffic (control) and were also used to assess pre treatment soil properties (Table 1). Soil bulk density was calculated as the mass of the oven-dried (105 °C over 72 h) soil sample divided by the volume of the cylinder. Soil bulk density was corrected for the >2 mm fraction, i.e. the volume and mass of >2 mm soil fragments were subtracted from the soil volume and dry mass, respectively. We found a significant and negative correlation between BD and WC (see Section 3.1 and Goutal et al. 2012) in the pre treatment sampling. Therefore, the day the forwarder drove on the soil, only initial WC was measured in the area where the soil was going to be disturbed by heavy traffic at 28 locations × 5 depths in CA (0–10, 10–20, 20–30, 30–40 and 40–60 cm) and at 32 locations × 3 depths (0–10, 10–30 and 30–50 cm) in AZ.

2.3.2 State following the forwarder traffic

Within 1 month after the forwarder traffic, two pits per block (three blocks per site) were dug in the trafficked area. In each pit, two soil samples per depth were collected to measure BD. We used 250 cm³ steel cylinders to collect undisturbed soil samples every 5 cm at 0–40 cm depth (with reference to the surface), i.e. in the silt loam layer, the soil layer the more likely to deform under the applied stress. The volume of the core sampler was expected to be enough to capture the entire soil porosity in the trafficked areas, where macro-porosity was assumed to be small. During sampling, we avoided the few very deep ruts accompanied with bulges at the edges that were formed by wheeling on the wettest soil conditions (plastic deformation; Williamson and Neilsen 2000; Ampoorter et al. 2010) and collected soil samples under the ruts averaging 5 cm depth.

2.4 Model simulation

2.4.1 Vertical stress distribution at the soil–tyre interface

One of the most important factors determining the performance of soil compaction models was found to be the upper boundary conditions, i.e. the soil–tyre contact area and the vertical stress distribution at this surface (Keller et al. 2007; Keller and Lamandé 2010). Consequently, we estimated the contact area and the distribution of normal stresses at the soil–tyre contact area using the model proposed by Keller (2005) that is incorporated in the SoilFlex model. We did all the simulations taking into account two load steps with four wheels (two front and two rear wheels), in order to mimic the stress applied in the field experiment.

2.4.2 Vertical stress distribution in the soil profile

The SoilFlex model is based on semi-analytical equations for stresses propagation in soil given by Söhne (1953). These equations include the so-called concentration factor that governs the decay pattern of stress with depth (for details, see Söhne 1953 and Keller et al. 2007). According to Défossez et al. (2003), we chose a concentration factor of 6, because our soils had a small density and were moist.

2.4.3 Soil stress–strain relationship

Calculations of soil deformation require soil mechanical parameters that characterise the stress–strain relationship. We chose the stress–strain relationship developed by O’Sullivan and Robertson (1996) as proposed in the SoilFlex model. The virgin compression line (VCL), and the recompression line (RCL) are given as (O’Sullivan and Robertson 1996):

$$ \mathrm{VCL}\,:\,v = N - {\lambda_n}\,\ln (p) $$

(1)

$$ \mathrm{RCL}\,:\,v = {v_{\mathrm{init}}} - \kappa\,\ln (p) $$

(2)

where: ν, soil-specific volume (dimensionless) = ρ _s/BD; ρ _s, density of solids (in grammes per cubic centimetre); p, mean normal stress (in kilopascals); N, specific volume (dimensionless) at p = 1 kPa; λ _n , compression index (in kilopascals); ν _init, initial specific volume calculated from the initial BD values (dimensionless); κ, recompression index (in kilopascals).

The mechanical parameters were measured by Saffih-Hdadi et al. (2009) using oedometer tests on remolded soil samples from the 10- to 30-cm soil layer of the AZ site. The measurements were not performed on intact soil cores, as the aim here was to calibrate soil mechanical parameters as a function of BD and WC, on the largest validity domain possible. The use of remolded samples is justified because we only use the parameters N and λ _n in our study. These two parameters define the VCL, which is measured at large stress levels well beyond the soil strength, and therefore, soil structure has only little impact (Keller et al. 2011). The samples were equilibrated at three WC (0.31, 0.25 and 0.2 g g⁻¹) corresponding to three matric potentials (−5, −33 and −100 kPa), and then samples at each initial water content WC_i were compressed to three initial BD values (1.1, 1.3 and 1.45 g cm⁻³) before undergoing oedometer tests (Saffih-Hdadi et al. 2009). For both sites, the mechanical parameters N and λ _n were estimated as functions of the initial measured values of BD and WC for each discriminated soil layer according to the data of Saffih-Hdadi et al. (2009). We used the same relationships for the CA than for the AZ site, because soil texture and SOC were similar (SOC, clay and silt content between 0.9–3.3, 16–36, and 40–64 % respectively for the AZ site and between 0.3–1.6, 11–23, and 59–74 %, respectively for the CA site). We chose to perform SoilFlex simulations on the medium and extreme initial conditions (Tables 2 and 3), i.e.:

• At the mean BD found before compaction associated with the mean WC during compaction of the soil layer considered (case A),

• At the largest BD value found before compaction associated with the smallest WC value during compaction (case B) and

• At the smallest BD value found before compaction associated with the largest WC value during compaction (case C).

Table 2 Initial soil parameters used for model simulations at the site of Azerailles

Table 3 Initial soil parameters used for model simulations at the site of Clermont en Argonne

Unfortunately, Saffih-Hdadi et al. (2009) did not measure the rebound and recompression index (κ) that is needed in the model of O’Sullivan and Robertson (1996). We estimated κ in two ways. (1) By using the pedo-transfer function of O’Sullivan et al. (1999) incorporated in the SoilFlex model, which is given as (O’Sullivan et al. 1999):

$$ \kappa = {\lambda_n}\left[ {0.119 - \left( {\frac{{0.082\ \mathrm{WC}}}{17 }} \right)} \right] $$

(3)

and (2) by calculating κ as:

$$ \kappa\,\frac{1}{3}\,{\lambda_n} $$

(4)

which was suggested by O’Sullivan and Robertson (1996) for wet soils. With these two approaches, we could evaluate the impact of κ on the simulated final (i.e. after wheeling) BD.

2.4.4 Model validation

Model validation was done by comparing simulated BD with measured BD. To be as close as possible to the field conditions when the driver of the forwarder created adjacent wheel tracks, we considered the simulated values across the width of the tyre, i.e. the 30 cm on either side of the wheel centre. For each 5 cm depth increment, we calculated the mean, minimum and maximum values of the predicted BD after wheeling across the wheel width.

The two replicates of BD measurements per depth and pit described in the previous section were taken at each side of the pit front (approximately one meter width). Therefore, the probabilities to have sampled at a given distance from the wheel track centre were considered as equal. By considering the SoilFlex predictions on this 60-cm wide stripe we assumed to have a modelled population as close as possible to the observed population.

2.5 Analysis of the effect of SOC on BD

To analyse the effect of SOC on BD of the undisturbed and trafficked soils, we sampled 10 randomly distributed soil samples per block × treatment × site from the 0- to 10-cm soil layer using 500 cm³ steel cylinders in April 2009. We sampled the 30 × 50-m areas of each treatment (undisturbed and trafficked) per block. After determining the field-moist sample weight to calculate WC, a sub-sample was taken to measure the BD as described in the previous section (2.3.1). The rest of the soil sample was air-dried and ground (sieve size of 250 μm) for determination of total carbon (C) and nitrogen (N) content using a CHN apparatus (ThermoQuest CE Instrument NCS 2500). The soil did not contain any carbonate; therefore, total carbon corresponds to the SOC content. The effect of SOC or C/N ratio, WC and treatment on BD was assessed through an analysis of covariance using the R software (version 2.11.1). The best model was selected after checking on the distribution of the residuals and by looking for the maximum coefficient of determination (R ²) and minimum residual standard error (RSE).

3 Results

3.1 Variability in initial BD and WC and observed soil compaction

The variability in BD and WC before heavy traffic was large (Tables 2 and 3) for both sites. For example, at 0–10 cm depth of each site, the maximal WC value was approximately twice as large as the minimal WC. This large variability in soil physical properties may have had an important impact on soil deformation following the forwarder traffic. For both sites, we found that BD was significantly correlated with WC at the sampling time (AZ: Pearson’s product–moment correlation of −0.47, t = −7.5, df = 194, p value = 2.3e−12; CA: Pearson’s product–moment correlation of −0.58, t = −8.8, df = 154, p value = 3.3e−15). Consequently, to take into account the effect of the variability in initial soil conditions on soil deformation, we considered the fact that WC and BD were not statistically independent. For example, the smallest value of BD before traffic was associated with the largest value of WC during traffic (Tables 2 and 3).

The soil of both sites was considered as highly sensitive to compaction, which was verified in situ. Despite the large variability in BD values, the effect of the forwarder traffic on BD was significant for the entire profile (0 to 35 cm depth) at both sites (Fig. 1).

Fig. 1

Soil bulk density measured before (control) and immediately after heavy traffic (trafficked); 18 and 12 replicates per depth in the Azerailles (AZ) site and 13 and 12 replicates per depth in the Clermont en Argonne (CA) site for control and trafficked, respectively. The grey-shaded areas and the bold line in the middle correspond to the range and to the mean of the measured bulk densities, respectively

3.2 Determination of SoilFlex _n and N parameters

Using the data of Saffih-Hdadi et al. (2009), we found the following relationships between N and _n , respectively, and BD and WC (coefficient of determination of 0.97):

$$ {\lambda_n} = 0.699 - 0.3575 \times \mathrm{BD} - 0.0029 \times \mathrm{WC} $$

(5)

$$ N = 6.324 - 2.249 \times \mathrm{BD} - 0.0296 \times \mathrm{WC} $$

(6)

Using these relationships, we took into account the variability in initial WC and BD in model simulations (Tables 2 and 3), and analysed its effect on soil deformation.

3.3 Stress simulation

The mean contact area modelled after the passage of four wheels (two front and two rear wheels) was 0.35 m² for both sites. The mean vertical contact stress was 156 and 153 kPa with a maximum value of 423 and 413 kPa under the edge of the tyres for the AZ and CA sites, respectively.

3.4 Effect of the rebound parameter on model estimations

The model overestimated BD for the 0- to 20-cm soil layer when using κ according to Eq. 3 (Fig. 2a). The mean of the simulated BD values for the 0- to 20-cm soil layer was of 1.67 and 1.43 g cm⁻³, with one of the observed values being of 1.39 and 1.30 g cm⁻³ for the CA and AZ sites, respectively. From 20 to 50 cm depth, the impact of the forwarder traffic on BD was slightly underestimated by the model. In the CA site, this tendency was larger: the differences between simulated and measured values were greater in comparison to the AZ site.

Fig. 2

Soil bulk density measured immediately following heavy traffic (field t ₀) and modelled considering mean initial conditions (case A, see text for details). a Simulations using Eq. 3; b simulations using Eq. 4. The grey-shaded areas and the bold line in the middle of the shaded area correspond to the range (min and max across the tyre width, i.e. 30 cm on each side of the wheel centre) and to the mean of bulk density values, respectively

The model-based estimations of BD were closer to the measured values for the surface layers of both sites when applying κ according to Eq. 4 (Fig. 2b). The improvement was particularly evident for the CA site. For the 0- to 5-cm soil layer, the difference between the mean of simulated and observed values decreased from 0.52 to 0.25 by using Eq. 4 instead of Eq. 3 for the CA site, whereas it decreased from 0.28 to 0.11 for the AZ site. For both sites, using Eq. 4 instead of Eq. 3 decreased the overestimation in the surface layer, but increased the underestimation in the deeper soil layers.

3.5 Effect of the variability in initial BD and WC on model-based estimations

Increasing initial WC and decreasing initial BD (case C) increased the range of BD values simulated after heavy traffic (Fig. 3). Yet, simulations never yielded as small values as we observed in the soil surface layer. For the AZ site, the simulations in case C predicted larger BD in the upper soil layers and smaller BD in the deep soil layers in comparison to case B (minimum initial WC and maximum initial BD). Bulk density predicted in case C using κ according to Eq. 3 was close to the observed values for depths below 20 cm.

Fig. 3

Soil bulk density measured immediately following heavy traffic (field t ₀) and modelled considering min and max initial conditions (cases B and C, see text for details). a Simulations using Eq. 3; b simulations using Eq. 4. The grey-shaded areas and the bold line in the middle of the shaded area correspond to the range (min and max across the tyre width, i.e. 30 cm on each side of the wheel centre) and to the mean of bulk density values, respectively

For the CA site, the simulations in case C yielded unreasonably large BD (Fig. 3a), and no deformation could be calculated using Eq. 4 in case C (Fig. 3b). For the CA site, soil mechanical parameters were not measured but estimated based on the relationships developed for AZ, as noted in Section 2.3.3. The initial WC was considerably larger for case C for CA than for any cases for AZ (Tables 2 and 3), which probably yielded inaccurate estimates of λ _n and N, and therefore inaccurate predictions of BD.

For the AZ site in case B and using Eq. 3, the model overestimated BD throughout the entire soil profile, whereas for the CA site the overestimation was restricted to the upper soil layers. For both sites in case B, using Eq. 4 improved model estimations in the upper soil layers, while the effect of κ on BD was marginal for depths below 20 cm.

3.6 Model performance

To assess the performance of SoilFlex in predicting soil compaction under forwarder traffic, we considered on the one hand the predictions from case A and on the other hand the mean of the simulated values from the three cases A, B and C. Based on the results presented in Fig. 2, we calculated soil deformation by using Eq. 4 for the 0- to 10-cm soil layer and by using Eq. 3 for the 10- to 50-cm soil layer (Fig. 4). That is, we used a larger rebound in the uppermost soil layer. Comparison of the predicted values to the observed ones yielded a root mean square error (RMSE) of 0.06 and 0.11 for the AZ and CA site, respectively. The lower performance of the model for the CA site was mainly due to the overestimated impact of the forwarder traffic in the surface soil layer even when using Eq. 4.

Fig. 4

Model-based predictions considering mean initial conditions (case A, see text for details), two passes of four wheels (two front and two rear wheels) and using Eq. 3 for the 0- to 10-cm soil layer and Eq. 4 for the 10- to 40-cm soil layer versus observed soil bulk density values before (control) and after heavy traffic (trafficked); 18 and 12 replicates per depth in the Azerailles (AZ) site and 13 and 12 replicates per depth in the Clermont en Argonne (CA) site for the control and trafficked areas respectively (means and standard deviation as error bars)

If we consider that the three simulation cases (cases A–C, respectively) were three locations on each site where we came to sample after heavy traffic, we would as measurement of the impact of forwarder traffic take the mean of the BD values per depth and site. Therefore, we calculated in the same way the mean of the three cases (AZ site) or of case A and B (CA site; unrealistic values for case C were discarded, as mentioned above). The RMSE was then 0.07 and 0.10 for the AZ and CA site respectively. These RMSE values are almost identical to the ones obtained for the case A simulation.

3.7 Effect of SOC on the impact of heavy traffic on soil BD

For both sites in the 0- to 10-cm soil layer of the control plots (C treatment), variations were best explained with a linear regression between BD and the logarithm (base e) of SOC and WC (Table 4).

Table 4 Results of the analysis of covariance accounting for effects of treatment (T; trafficked), soil gravimetric water content (WC; in g 100 g⁻¹), and soil organic carbon content (SOC; in g 100 g⁻¹) on soil bulk density (in grams per cubic meter)

Bulk density increased with decreasing SOC (Fig. 5a) and with decreasing WC (Table 4). No effect of the C/N ratio could be observed. Soil water content increased with increasing values of SOC (Table 5; Fig. 5b).

Fig. 5

Partial relationships between soil bulk density and ln(SOC) while the second covariate is held constant at its arithmetic mean (a) and relationships between gravimetric water content (WC) and ln(SOC) (b). For each graph, C correspond to the control treatment (undisturbed areas) and T to the trafficked treatment; 30 samples per site and treatment. Sampling in April 2009, 1 (CA site) to 2 (AZ site) years after heavy traffic

Table 5 Results of the analysis of covariance accounting for effects of treatment (T; trafficked) and soil organic carbon content (SOC; in g 100 g⁻¹) on soil gravimetric water content (WC; in g 100 g⁻¹)

The effect of SOC on the increase in BD after wheeling was expected to be large because of its effect on initial BD and WC. According to Eqs. 5 and 6, λ _n and N are affected by both BD and WC, but the relationship between SOC and mechanical properties is complex: while SOC decreases BD, it increases WC. However, the effect of BD on λ _n (Eq. 5) and N (Eq. 6), respectively, is about one order of magnitude larger than the effect of WC on λ _n (Eq. 5) and N (Eq. 6), respectively. Therefore, both λ _n and N increase with increasing SOC. This results in a decrease of compaction at small stresses, but in an increase of soil compaction at large stresses such as during loading with the forwarder. However, no significant correlation between the effects of SOC and treatment could be detected, i.e. the increase in BD from control to trafficked soil did not increase with increasing SOC (Table 4). This is also seen from the fact that the linear regression lines between BD and log SOC for the different treatments were parallel (Fig. 5a). Hence, the larger λ _n (larger compressibility) at larger SOC seems to be compensated with a larger κ (larger rebound).

4 Discussion

4.1 Rebound may be increased in surface layers of forest in comparison to agricultural soils

For both sites and for the three initial soil conditions, no model simulation could predict the smallest BD found in the field at 0–5 cm depth after the forwarder traffic and only a slight increase in BD in the 0- to 5-cm soil layer in comparison to the increase at 5–10 cm depth. One explanation could be that all the analytical procedures used in the SoilFlex model were developed for agricultural soils and do not consider the buffering effect of the forest floor and the change in soil mechanical properties due to the large organic matter content and persistent root mat. The presence of organic matter may have increased the rebound (Soane 1990). Indeed, a clear improvement of predictions was observed for the surface soil layers when increasing the rebound by using Eq. 4 instead of Eq. 3. Besides, for both sites we found, in the 0- to 10-cm undisturbed soil layer, a positive correlation between SOC and WC, and a negative correlation between SOC and BD (Fig. 5). Goutal et al. (2012) also found at AZ and CA a relationship between BD and WC which was related to soil swelling and to the fact that WC at sampling is well correlated to soil constituents. These results support the fact that larger SOC is associated with larger WC but also with smaller BD (swelling due to large WC and effect of SOC on BD). According to the effect of BD and WC on soil susceptibility to compaction (Saffih-Hdadi et al. 2009), a greater amount of SOC should have increased soil compaction during forwarder traffic. Yet, the statistical evaluation of the relationship between BD and SOC did not show a dependence of the treatment effect on SOC. At AZ, the linear regression lines between BD and ln(SOC) of the undisturbed and the trafficked soils are shifted but parallel to each other (Fig. 5a). The impact of heavy traffic (i.e. the shift) was constant whatever the SOC. Williamson and Neilsen (2000) found similar results between BD and organic matter before and after 15 passes of a laden logging machine. Schäffer et al. (2008) also found similar results concerning the respective effects of SOC and traffic on textural and structural porosities of agricultural soils. At CA, the same conclusion could be drawn, except that the impact of heavy traffic on BD (shift) was not significant (Fig. 5a). We only measured BD and were not able to analyse the effect of SOC on textural and structural pore spaces (Tranter et al. 2007; Pereira et al. 2007; Boivin et al. 2009; Schäffer et al. 2008). For example, after accounting for the SOC content effect on BD, the relationship between WC and BD was affected by treatment for the AZ site (the slope between BD and WC was −0.2 and −1.4 for the trafficked and control plots, respectively). This may be related to a change in swelling-shrinking properties after compaction (Goutal et al. 2012) and to a modification of the soil structure (Boivin et al. 2009; Schäffer et al. 2008). At CA, the relationship between SOC and WC was affected by treatment, which could originate from different impacts on the textural and structural pore space. The effect of SOC on the resistance to stresses may have been different for the texture- and structure-related pore spaces, resulting in similar BD after heavy traffic for large SOC and small SOC locations. According to the results of Schäffer et al. (2008), this assumption can be rejected as they found no effect of SOC on the impact that heavy traffic had on soil textural and structural porosities. These results underlined an additional effect of SOC on soil mechanical behaviour besides its effect on initial BD and WC. It was namely possible that soil with large SOC content was highly compressible but also recovered more after unloading (rebound) (Kuan et al. 2007). That is, SOC increases compressibility (i.e. λ _n ; Eq. 1) as well as resilience (i.e. κ; Eq. 2). This seems to agree with our simulations, where the introduction of an increased rebound improved the model performance. Therefore, a future improvement of analytical soil compaction models could be to take into account the effect of SOC on soil mechanical parameters, especially on the rebound parameter, but also on stress propagation.

In general, SoilFlex predictions tended to over-estimate the impact in the upper soil layers and to underestimate it in the deeper soil layers. Therefore, even if forest floor litter and/or soil organic matter decreased the impact of the forwarder on the soil surface layers, stresses were transmitted to the deeper soil layers. In our study, the 0- to 10-cm soil layer had a very small density in comparison to the rest of the soil profile; the current one-layer analytical approach may therefore be a possible cause of model accuracy loss (Keller and Lamandé 2010). Furthermore, WC of the soil surface layers were at or slightly beyond the wet end of the range at which mechanical properties were measured, and therefore, estimation of mechanical properties (Eqs. 5 and 6) may have been less accurate for these layers as compared with the deeper soil layers. According to the values of λ _n calculated thanks to Eq. 5, the pédo-transfer function of O’Sullivan et al. (1999) (Eq. 3) yields negative values for WC > 25 g 100 g⁻¹, another possible reason why model simulations with Eq. 3 do not work well in the surface layers. Besides, because of the wet soil surface layer, it is possible that a combination of compaction and shear (with potential loosening) occurred at 0–5 cm depth, conditions which are not well simulated by SoilFlex (or any other soil compaction model). It may also be possible that predictions of stresses (and therefore BD) were less accurate and included an uncertainty that is greater close to the surface and decreases with depth, because the contact stresses were estimated and not measured. Besides, the prediction equations by Keller (2005) were developed from measurements using agricultural tyres on arable soil and the small density due to large SOC at soil surface may have also resulted in an increased soil–wheel contact area and thus smaller contact stresses. We could not test this hypothesis using SoilFlex as the model developed by Keller (2005) to describe stress distribution at the soil–wheel interface did not account for soil properties.

4.2 Evaluation of model performance

The predicted mean vertical contact stress fell into the range of values measured by Horn et al. (2007) for forest machinery. Taking into account the spatial variability of initial soil conditions did not improve model performance. Model simulations performed better for the AZ than for the CA site, because the relationships between soil mechanical and physical properties were measured for AZ but only estimated for CA. Different authors (Défossez et al. 2003; Keller et al. 2007) also found that the simulations of soil deformation under wheeled traffic are very sensitive to the assessment of soil mechanical properties. Considering the minimum BD and maximum WC at the time of traffic yielded unrealistic results for the surface layer of the AZ site when using Eq. 3, and for all depths for the CA site. The mechanical parameters N and λ _n calculated were probably not accurate in these cases. Indeed, Saffih-Hdadi et al. (2009) did not measure mechanical parameters beyond a maximum WC of 31 %.

The soil was wetter at the time the forwarder drove on the CA site in comparison to the AZ site. The greater improvement of predictions by using the rebound parameter estimation for wet soils (Eq. 4) for the CA site in comparison to the AZ site could be due to these very wet conditions. Yet, even by using Eq. 4, the model over-estimated the impact of the forwarder on the soil surface layer of the CA site, unlike the AZ site. AZ and CA depicted nearly identical SOC contents in soil surface layers (2.4–4.7 and 2.2–4.0 at 0–10 cm depth for the AZ and CA sites, respectively), but as they did not have the same WC at the time the forwarder drove on the site, the effect of SOC may have differed. Therefore, the very small impact of heavy traffic on the soil surface layers of the CA site in comparison to the theoretical expectations may be related to the specific role of SOC in wet conditions (Soane 1990; Mosaddeghi et al. 2000; Pereira et al. 2007). Another explanation could be that the 0–10 cm soil layer was very wet at the time of traffic and could flow (Williamson and Neilsen 2000; Ampoorter et al. 2010), whereas the minimum 0.1 g g⁻¹ drier 10–50 cm soil layer did not. It resulted in a very small impact on BD in the surface layer with only shallow ruts (minority of ruts exceeding 5 cm depth), in accordance with Mosaddeghi et al. (2000).

5 Conclusions

For the AZ site, we achieved the best model performance (RMSE = 0.06) considering the mean initial soil conditions (mean initial BD and WC) and the rebound parameter (κ) estimation for wet soil (O’Sullivan and Robertson 1996) from 0 to 10 cm depth and the estimation of O’Sullivan et al. (1999) from 10 to 50 cm depth. The model estimations predicted well the range of BD values found after forwarder traffic. Therefore, for soil conditions (BD and WC) in the range where soil mechanical parameters were measured, SoilFlex yielded accurate predictions even for forest soils with a non-negligible gravel and root content. Taking into account the variability of soil initial conditions did not improve model performance.

For the CA site, the wetter soil at the time of the compaction and the lack of accurate assessment of mechanical properties decreased the quality of the predictions. Yet using the O’Sullivan and Robertson (1996) estimate for κ from 0 to 10 cm and the O’Sullivan et al. (1999) one from 10 to 50 cm while considering the mean initial conditions led to satisfactory predictions (RMSE = 0.11).

For both sites, large SOC content was associated with large WC and small BD, yet the change in BD caused by heavy traffic was not found to be influenced by SOC at 0–10 cm depth. Therefore, soil displaying larger SOC content may have been more compressible due to larger WC and smaller BD, but probably recovered more than soils poor in SOC. Indeed, the soil deformation predictions for the 0- to 10-cm soil layer of both sites were initially poor, but improved when we increased the rebound parameter. Analytical models of soil compaction need further development to take into account the (mechanical) properties, soil heterogeneity and structure of the litter and organic layers of forest soils.