Characterization and prediction of hydraulic properties of traffic-compacted forest soils based on soil information and traffic treatments

Abstract

Key message

The hydraulic properties of compacted and rutted soils were evaluated through in-situ infiltration experiments and predicted based on soil texture class and traffic treatments. A significant decrease in saturated soil water content and soil hydraulic conductivity at saturation was observed. The resulting soil hydraulic parameters, when integrated into a soil water transfer model, effectively simulated water dynamics in these impacted forest soils, providing a crucial first step toward developing decision support tools for real-time trafficability. This approach can assist forest managers in minimizing the extent of soil compaction.

Context

To overcome trafficability issues of forest soils induced by heavy logging machinery, planning support tools are needed to determine suitable soil moisture conditions for traffic.

Aims

This study aimed to identify the soil properties that differ significantly between undisturbed and compacted soils and to provide several estimation tools to predict the hydraulic properties of compacted soils beneath the skid trails.

Methods

Four hundred seventeen water infiltration tests were conducted on 19 forest sites, mostly in North-eastern France, and analysed with the BEST method to estimate the hydraulic properties of the skid trails and undisturbed soils. The hydraulic properties of the skid trails were predicted thanks to linear mixed effect models using a bulk treatment effect, a site effect, or a skid trail degradation score. The predicted hydraulic properties were tested using a water flow model to assess their relevance regarding the prediction of water dynamics in skid trails.

Results

The compaction effect was only significant for the logarithm of the hydraulic conductivity at saturation (log₁₀(K_sat)) and the soil water content at saturation (θ_sat). For the skid trails, θ_sat was reduced by - 0.02 and − 0.11 m³m⁻³ in the 0 − 10 cm and 15 − 25 cm layers respectively, compared to undisturbed topsoil (0 − 10 cm). log₁₀(K_sat) was reduced by − 0.38 and − 0.85 for skid trails in the 0 − 10 and 15 − 25 cm soil layers respectively, compared to undisturbed topsoil. The use of a pedotransfer function, in replacement of water infiltration tests, and their combination with the same correction coefficients proved to efficiently simulate the difference in water dynamics between skid trails and undisturbed forest soils.

Conclusion

Estimation of soil hydraulic properties based on in situ water infiltration experiments proved efficient to simulate water dynamics in compacted and rutted forest soils. Yet, further studies are needed to identify the most adapted pedotransfer function to forest soils and to test the generalisation of our findings in different conditions, especially deeply rutted soils (rut depths > 12 cm).

1 Introduction

Soil compaction induced by human activity is recognized as one of the major threats for soil quality. Forest soil compaction has increased with heavy logging machinery like harvesters, skidders and forwarders (Greacen and Sands 1980). The consequences of soil compaction are now well known, as changes in soil structure with the increase of dry bulk density (Kozlowski 1999; McNabb et al. 2001; Ampoorter et al. 2007) and the reduction of soil porosity, especially macroporosity (McNabb et al. 2001; Frey et al. 2009). The loss of large pores considerably affects capillarity and water retention, soil aeration (Ballard 2000; Frey et al. 2009), and hydraulic conductivity and infiltration capacity (Kozlowski 1999). Reduction of soil water content and hydraulic conductivity at saturation were observed (Kozlowski 1999; Williamson and Neilsen 2000; Schack-Kirchner et al. 2007), with at the same time, increased water contents under dry conditions due to increased capillarity. Therefore, compacted soils remain wet over longer periods than undisturbed soil, leading to a reduction in the period available for machinery traffic, since the soil-bearing capacity is strongly reduced for wet soils (McNabb et al. 2001). For forest management, this leads to strong constraints on harvesting logistics.

Soil compaction has a huge and long-lasting impact on forests, especially it leads to mean decreases in total biomass of about 20% and in root depth of about 30% (Mariotti et al. 2020), with a recovery taking at least 50 years (Mohieddinne et al. 2019). Therefore, preventing the occurrence of soil compaction is of paramount importance for forest managers. Forest managers can act at two levels. The first induces the reduction of the area dedicated to forestry engine circulation with the creation of specific skid trails. Forest machineries must drive on these skid trails leaving the remainder of the plot free from heavy machinery circulation and thus minimizing the forest area impacted by soil compaction. The second operates when soil-bearing capacity can support heavy forest machineries, which limits soil compaction and maintain long-term trafficability of the skid trails. Occurrence of unfavourable traffic conditions under skid trails can be increased by water accumulation in the ruts and slower soil drying due to compaction. Compaction and rutting effects on water dynamics are difficult to account for since they depend on soil type and soil degradation intensity. So far, there are only few characterization support-tools of soil degradation mostly documented by qualitative visuals indicators like rut depth classification (Heninger et al. 2002; Frey et al. 2009) and only a few soil bearing capacity prediction support-tools (Ziesak 2003; Goutal et al. 2013; Salmivaara et al. 2021; Vega-Nieva et al. 2009). Goutal et al. (2013) have modelled the mechanical properties of two forest soils to predict their deformation with water content as input. Ziesak (2003)’s support-tool also requires soil moisture measurements by the forest manager, which is difficult to put into practice. Vega-Nieva et al. (2009) and Salmivaara et al. (2021)’s support-tool use water flow modelling to predict soil moisture in real-time data then improving accessibility to forest managers. Characterizing soil hydraulic properties, which are inputs of soil water transfer models, is then a crucial issue, especially in already compacted soils as found under skid trails.

Soil hydraulic properties are rarely characterized by forest practitioners, especially under skid trails. As a matter of fact, measurements are time-consuming and soil heterogeneity requires large sampling designs. Alternatives are to use either pedotransfer functions (PTF) based on soil proxies like soil texture or soil organic matter content for undisturbed soils (Wösten et al. 1999; Schaap et al. 2001; Tóth et al. 2015), bulk density for compacted soils (Assouline 2006) or visual proxies as rut depth classification for skid trails (Heninger et al. 2002; Frey et al. 2009). PTF are mainly developed to estimate soil hydraulic properties for undisturbed soils but may not be suitable for compacted forest soils, especially for deeply rutted soils displaying small increase in bulk density and large decrease in hydraulic conductivity (Frey et al. 2009). Some PTFs use bulk density variation, which is an indicator of compaction, as a predictor of soil hydraulic properties (Wösten et al. 2001; Schaap et al. 2001; Assouline 2006; Tóth et al. 2015). However, databases used to develop these PTFs contain only a small fraction of forest soils, that are characterized by high organic matter content (Wösten et al. 1999; Weynants et al. 2013; Al Majou et al. 2008) and a very small proportion of compacted soils. Besides, they are not considering structural changes induced by rutting and related impacts on the soil hydraulic properties while Pachepsky and Rawls (2003) claimed that the use of information on soil structure would improve PTF performance.

The overall objectives of this study were to provide a reference for soil hydraulic properties of compacted soils under skid trails and design tools to estimate such properties from proposed proxies.

Our hypotheses were that (i) soil hydraulic properties are impacted by compaction and rutting, as the result of the effect of compaction on pore size distribution, (ii) visual indicators and proxies may help to predict such impact, the higher the visual score of skid trail degradation the lesser the hydraulic conductivity, and (iii) such indicators and proxies can be used to predict the hydraulic properties from undisturbed soils and use them for modelling water dynamics in compacted and rutted forest soils.

To reach our objectives, we firstly determined the hydraulic properties of numerous forest soils using the Beerkan method (single-ring water infiltration experiments) and treating the data with the Beerkan Estimation of Soil Transfer parameters (BEST) method (Braud et al. 2005; Lassabatere et al. 2006). Then, we compared the data obtained under skid trails with those obtained in close undisturbed soils. Secondly, we relied on these data to compare two approaches to predict the soil hydraulic properties of skid trails based either on water infiltration measurements or using PTF estimations. We tested, for each approach (either BEST or PTF ), the relevance of including a visual score of skid trail degradation intensity to estimate soil hydraulic properties of compacted soils. Lastly, the predicted hydraulic properties were used to model water flow through a profile of soil and compared to observed data.

2 Material and methods

2.1 Study sites

Nineteen plots located in 16 forests were studied (Table 1, Fig. 1). Firstly, we selected 12 plots located in Northeastern France affected by soil trafficability issues because of rainfall occurrence during forest exploitation periods and poorly draining soils. Each selected site satisfied the following criteria: plain deciduous stand, fine textured soils sensitive to compaction and presence of skid trails already trafficked. Then, the initial dataset was completed by four additional plots located in Central and Southern France to increase the diversity of soil texture and soil degradation level in the database. Time since the last heavy traffic was not known. Then, the selection of experimental sites was based on time before a forthcoming exploitation considering short times before the next harvesting (1 to 2 years) in order to avoid any need for detailed records on former traffic events and natural recovery processes. The experiments were carried out on both undisturbed soils hereafter named Control (C) treatment, and circulated trails hereafter named Trafficked (T) treatment. Traffic intensity in terms of used machinery, number of machines passes and weather conditions at the time of traffic was not known for 14 plots. The 2 plots where traffic intensity was known and soil recovery was monitored are the two experimental sites described by Goutal et al. (2012), and named hereafter “AZ25” and “CA”. The timeframe between two successive forest operations are about 7–15 years for the studied forests; therefore, the measurements at AZ25 and CA were done in the same timeframe.

Table 1 Name, location, stand composition of the studied plots and intensity of degradation of the skid trails

Fig. 1

Location of studied plots

2.2 Characterization of skid trail degradation

In order to remain in an operational context, characterization of the intensity of the skid trail degradation was only performed on the basis of qualitative visual indicators. No reference was made to the traffic records or details on machinery used during previous forestry works. Indeed, these data are rarely available for forest operators. We adapted the soil disturbance classification developed by Heninger et al. (2002). This operational classification was based on several observations easily accessible with an auger and a spade (Fig. 2):

Fig. 2

Field soil disturbance classification grid. Rut depth is measured below the original soil surface and corresponds to the deepest soil depression in the skid trail at the measurement’s location. For each criterion the score is determined according to rut depth, soil structure and the apparition depth of redoximorphic features. The sum of scores yields the degradation class of the skid trail: low, moderate and severe degradation class correspond to a sum of scores between 1 and 3, between 4 and 6 and between 7 and 11, respectively

•The soil structure at 0–10 and 15–25 cm observed in the C and T treatments.

•The depth of redoximorphic features indicating poor aeration conditions and water ponding. The difference between the depths of redoximorphic features measured on C and T treatments was calculated and split into 4 categories (0,]0–10 cm],]10–30 cm], ≥ 30 cm).

•The rut depth measured on the T plot indicating soil compaction and/or soil creep, split into 3 categories ([0–5 cm],]5–10 cm], ≥ 10 cm).

A score of disturbance between 0 and 3 was calculated for every criterion and then summed to get a global indicator. For each criterion, 0 was noted when the criterion in the T-treatment was equal to that of the C-treatment. When a difference could be observed between the T-treatment and the C-treatment, the score was comprised between 1 and 3, 1 being the lowest and 3 being the highest level of disturbance. A skid trail was considered as slightly disturbed if the sum of the scores was between 1 and 3, as moderately disturbed between 4 and 6, and as severely disturbed between 7 and 11. These three categories involved respectively 6, 23 and 12 forest sites.

2.3 Characterization of soil hydraulic properties

Four to eight Beerkan experiments were performed per plot, per depth, i.e. soil layer (0–10 cm and 15–25 cm), and per treatment. In order to have at least three successful Beerkan experiments per situation, their number was increased on some plots because some of the tests lead to very small infiltration or even zero infiltration due to extremely low hydraulic conductivity, especially in the skid trails. In these cases, the experiments were aborted and considered as failed, and additional experiments were performed nearby. For each Beerkan experiment, a soil core sample was extracted to determine the bulk density (ρ_b) and the initial water content θ₀. A composite sample was collected for each unique combination of site, depth and treatment in order to analyse the particle size distribution, pH and organic matter content (Table 2).

Table 2 Mean soil properties: clay and silt contents, USDA (US Department of Agriculture) classification system (Hartshorne and Dicken 1935) and HYPRES (European classification system) (Wösten et al. 1999) textural classes, soil organic carbon content (SOC), pH in water (1:5 in volume)

In situ Beerkan infiltration tests consist in recording the infiltration time of a succession of fixed volumes of water poured into a simple annular ring (Braud et al. 2005; Angulo-Jaramillo et al. 2016). According to the measurement campaigns, the radius of the cylinders varied between 50 and 98 mm while the water volumes poured into the ring varied between 60 and 75 cm³ corresponding to a water height between 2.5 and 8 mm (Martin et al. 2024). Note that BEST methods account for the difference in experimental conditions (ring size, initial water content, etc.) so that the estimated hydraulic parameters are independent of experimental conditions as long as the measurements properly describe the transient and steady states of water infiltration.

Setting aside the aborted Beerkan experiments, five alternative BEST methods were used to treat the 399 successful infiltration experiments. The unsaturated hydraulic properties of soils involve the two following functions: the unsaturated hydraulic conductivity curve K(θ) or K(h) (m s⁻¹) representing the soil capacity to conduct water as a function of the soil volumetric water content (VWC) θ (m³ m⁻³), or the soil matric potential h (m), and the water retention curve h(θ) representing the soil capacity to hold water by capillarity. These soil water retention and hydraulic conductivity curves were estimated with the BEST method using the van Genuchten relationship (van Genuchten 1980) with the Burdine condition (Burdine 1953) to define h(θ) (Eq. 1) and the Brooks and Corey relationship (Brooks and Corey 1964) to define K(θ) (Eq. 2).

$$S_e=\frac{\theta-\theta_r}{\theta_{sat}-\theta_r}=\left(1+\left(\frac h{h_g}\right)^n\right)^{-m}$$

(1)

$$ \frac{K\left(\theta\right)}{K_{sat}}=\left(\frac{\theta-\theta_r}{\theta_{sat}-\theta_r}\right)^\eta\text{with}\;\eta=\frac2{n m}+2+p$$

(2)

where θ_sat and θ_r are respectively the soil VWC at saturation and the residual VWC, h_g is the scale parameter for water pressure head (m, < 0), K_sat is the hydraulic conductivity at saturation (m s⁻¹), p is a tortuosity parameter fixed at p = 1 for Burdine’s condition, and n, m, and η are shape parameters. BEST methods estimate the shape parameters from the particle size distribution using the specific pedotransfer functions described in Lassabatere et al. (2006). θ_sat is deduced from the bulk density ρ_b measurement and the particle density (ρ_s = 2.65 g cm⁻³): θ_sat = 1 − ρ_b / ρ_s and θ_r is assumed to be null (Lassabatere et al. 2006).

The original BEST method was developed by Lassabatere et al. (2006) and then improved by Yilmaz et al. (2010), Bagarello et al. (2014), Di Prima et al. (2021), and Yilmaz et al. (2022). It is based on quasi-exact implicit equation developed by Haverkamp et al. (1990) for 1D water infiltration and extended to 3D for water infiltration into single rings (Smettem et al. 1994). Haverkamp et al. (1994) proposed explicit approximate expansions for the description of the transient and the steady states, with quite accurate precisions (Haverkamp et al. 1994; Lassabatere et al. 2009). BEST methods are easier to implement, less time-consuming than laboratory measurements and thus more replicates can be performed. They were successfully applied in many operational contexts (Angulo-Jaramillo et al. 2019). In this study, the cumulative infiltration curve I(t) obtained from Beerkan experiment was analysed with each of the five BEST methods to estimate the hydraulic parameters: BEST-slope (Lassabatere et al. 2006), BEST-intercept (Yilmaz et al. 2010), BEST-steady (Bagarello et al. 2014), BEST-WR (Di Prima et al. 2021) and BEST-WR-3 T (Yilmaz et al. 2022). Note that these methods provide different estimates for only the two scale parameters K_sat and h_g. This last parameter, h_g, is estimated from the previous estimations of the saturated hydraulic conductivity, K_sat, and the sorptivity, S (Lassabatere et al. 2006).

The first three methods are based on the same mathematical framework, but their fitting process differ (Martin et al. 2024). The two last methods, BEST-WR (Di Prima et al. 2021) and BEST-WR-3 T (Yilmaz et al. 2022), have the same fitting process as BEST-slope but adapt the mathematical framework to water repellency. The BEST-WR and BEST-WR-3 T methods generated lower fitting errors for water-repellent soils but also for hydrophilic or partially hydrophobic soils (Martin et al. 2024) as noticed by Di Prima et al. (2021) and Yilmaz et al. (2022).

We calculated the relative error and the bias to assess the quality of the fit for each infiltration experiment and for the 5 methods (Martin et al. 2024). The goodness of fit was based on the adequacy of the general shape of the model and the relative error as suggested by Lassabatere et al. (2006). All the fits with relative errors smaller than or equal to 5% were automatically kept (Lassabatere et al. 2006). All the fits with relative errors higher than 10% were automatically discarded. Between these two thresholds, some fits displayed quite accurate shapes without systematic over- or underfitting. Therefore, we chose to remove infiltration experiments showing a strong bias while keeping the others. The FTG and BTG sites were excluded from the dataset because no infiltration experiments could be selected in the T-treatment on the basis of the criteria defined. Then, the dataset was summarised per Beerkan experiments by calculating the arithmetic mean of the hydraulic properties estimated from the all the selected BEST methods as suggested by several authors (Lassabatere et al. 2019). All details on the fitting process can be found in the companion data paper (Martin et al. 2024).

2.4 Integration of aborted experiments into the dataset

Eighteen infiltration experiments over the 417 (4.3%) performed in our study were aborted because of very small infiltration rates, i.e. requiring more than 40 min for the infiltration of the first volume of water or more than 4 h for the infiltration of the five first volumes without attainment of the steady state. All experiments performed on C-treatment were successful while 8.4% of experiments performed on T-treatment were aborted. The aborted tests are equally distributed on the two layers (0–10 cm or 15–25 cm) and mostly located on the AZ, POU2 and POU3 sites. We tried to explain these aborted experiments with soil properties as initial water content, particle size distribution, organic matter, bulk density and soil disturbance class but statistical results were not conclusive (Martin et al. 2024).

In the case of aborted experiments, the BEST methods could not be used because of too few points to describe both the transient and steady states of the water infiltration. Therefore, parameters of the soil hydraulic curves could not be estimated for these runs corresponding to very low values of saturated hydraulic conductivity. Simply removing these aborted experiments would have led to an underestimation of the compaction effect on K_sat, since only the T treatment was concerned by experiment abortion, i.e. very low permeability resulting from soil compaction. To overcome this potential statistical bias, we chose to develop an empirical method to estimate K_sat for these aborted experiments.

We first calculated a K_sat majorant defined as the minimum value of K_sat allowing for infiltration of the first water volume in 40 min. Such a criterion was defined arbitrary but corresponds also to practical considerations on the field. For that purpose, cumulative infiltrations were simulated with the Haverkamp model (Haverkamp et al. 1994) with different values of K_sat. To set up this model, we used (i) the field conditions of the aborted infiltration experiments to set initial water content and bulk density, and (ii) the averaged values of sorptivity calculated from successful infiltration experiments (Fig. 3). In Fig. 3, the horizontal line corresponds to the first water height poured into the cylinder and the vertical line corresponds to the maximum infiltration time (i.e. 40 min = 2400 s). By doing so, we assumed that the values of sorptivity of these soils were not drastically different from the other soils. With K_sat > 5 10⁻⁴ mm s⁻¹, all experiments should have led to success whereas when K_sat < 10⁻⁴ mm s⁻¹ all experiments should have failed. This range was therefore considered relevant for the definition of the majorant of K_sat for aborted experiments.

Fig. 3

Cumulated infiltration curves of aborted infiltration experiments simulated with the Haverkamp equation (Haverkamp et al. 1994). The horizontal line corresponds to the first water height poured into the cylinder and the vertical line corresponds to the maximum waiting time for the first poured water volume. Intersection of the two lines gives an idea of the value of K_sat on skid trails where the infiltration experiments were aborted

Then, we built four different datasets in order to test the reliability of this empirical method to decrease analysis bias with regard to low-permeability soils. Dataset 1 was built with only results obtained from non-aborted experiments: this dataset is expected to produce the worst statistical bias. Datasets 2, 3 and 4 were built with results obtained from non-aborted experiments and a fixed value of K_sat for all aborted experiments. Values of K_sat were always lower than the majorant with fixed values set at 10⁻⁴ mm s⁻¹ for dataset 2, 10⁻⁵ mm s⁻¹ for dataset 3, and 10⁻⁶ mm s⁻¹ for dataset 4.

2.5 Statistical analysis

In order to provide (i.) the most reliable reference possible on the hydraulic properties of compacted forest soils (skid trails) and (ii.) estimation tools for soil hydraulic characterization, we characterized the effects of several in situ measured factors (treatment, depth, and intensity of skid trail degradation) on the hydraulic parameters defining the water retention and the hydraulic conductivity curves. For this purpose, we used different statistical models. Prior to statistical analysis, the hydraulic conductivity at saturation and the scale parameter h_g were log transformed because of their skewed distribution and the heteroscedasticity of their related residuals (i.e. contrasting standard deviations between groups).

The SITE_mod model considered the effects of two fixed effects “treatment” and “depth” and their interaction “treatment:depth”, while considering that experimental sites could produce random variations of the treatment and depth effects. The random effect represents the spatial variability between undisturbed sites (R_s in Eq. 3) and spatial variability of the intensity of the soil compaction during the traffic, as the consequence of contrasting soil sensitivity to compaction and different traffic conditions (R_sj in Eq. 3). The “lme4” R package (Bates et al. 2015) was used to compute the different terms of the SITE_mod model given in Eq. 3.

$${\widehat y}_{ijk}=\mu+\alpha_i+\beta_j+\alpha_i:\beta_j+R_s+R_{sj}+\varepsilon_{ijk}$$

(3)

where ${\widehat{y}}_{ijk}$ is the predicted value for sample k at the ith depth (0–10 cm or 15–25 cm) and jth treatment (Control or Trafficked), µ the intercept corresponding to the C-treatment at 0–10 cm, ${\alpha }_{i}$ and ${\beta }_{j}$ are the fixed effects with ${\alpha }_{i}$ the depth effect (${\alpha }_{i}$= 0 at 0–10 cm) and ${\beta }_{j}$ the compaction effect (${\beta }_{j}$= 0 for C-treatment), ${\alpha }_{i}:{\beta }_{j}$ is the interaction effect (traffic effect at 15–25 cm depth if significantly different from $\mu +{\alpha }_{i}+{\beta }_{j}$), ${R}_{s}$ is the site random effect, ${R}_{sj}$ is the site random effect for treatment j and ${\varepsilon }_{ijk}$ are the residuals.

In order to increase model reliability, we tested a model where the treatment effect was replaced by the soil disturbance class (i.e. visual indicator of skid trail degradation), named hereafter as PROX_mod model.

$${\widehat y}_{ilk}=\mu+\alpha_i+\gamma_l+\alpha_i:\gamma_l+R_s+\varepsilon_{ilk}$$

(4)

where ${\widehat{y}}_{ilk}$ is the predicted value for sample k at the ith depth (0–10 cm or 15–25 cm) and l level of the visual indicator of skid trail degradation, µ the intercept corresponding to the control-treatment at 0–10 cm, ${\alpha }_{i}$ is the depth effect (${\alpha }_{i}$= 0 at 0–10 cm), γ_l is the coefficient for the level l of visual indicator of skid trail degradation (γ_l = 0 for the C-treatment), ${\alpha }_{i}:{\gamma }_{l}$ is the interaction effect between depth and visual indicator of skid trail degradation, R_s is the site random effect (related to the between-site variability of soil properties for the C-treatment at 0–10 cm) and ${\varepsilon }_{ilk}$ are the residuals.

Each model was simplified by keeping the most significant factors according to the likelihood ratio chi-squared test applied for successive nested models. For each retained model, the bootstrap method (Efron and Tibshirani 1993) was used to calculate mean and standard deviation of the model coefficients. Sampling with replacement was done 1000 times, randomly among the data set, each random sample having the size of the initial data set. For each random sample, the retained model was fitted. Model coefficients calculated on each random sample were then summarised by the average and the standard deviation.

2.6 Comparison of the different approaches used to predict soil hydraulic parameters of skid trails

Development of SITE_mod and PROX_mod models aims to simplify the estimation of hydraulic properties of compacted soils, which are required for the simulation and prediction of water flow. If the estimation of soil hydraulic parameters is not perfect, one might wonder if the impact of compaction on soil water dynamic may be properly modelled on the basis of our results. To answer this question, we used a mechanistic model of water flow (Lafolie 1991; Findeling et al. 2007) to simulate the water dynamics in compacted soils and in undisturbed soils. Modelling of the soil water dynamics was performed with the Virtual Soil (VSoil) modelling platform developed by the National Institute of Agricultural Research (https://vsoil.hub.inrae.fr/). Our model, named forest_soil_practicability_prevision_tool, is based on the PASTIS model (Prediction of Agricultural Solute Transfer In Soils) developed by Lafolie (1991) and improved by Findeling et al. (2007) and Martin (2019). It is a one-dimensional mechanistic model of soil–plant-atmosphere system able to simulate the water and heat transfer in soils. It was proved to be robust and reliable to simulate soil water transfer (Findeling et al. 2007) on a physical basis using Richards equation (Richards 1931) and Fourier’s law to represent soil water and heat flows, respectively.

Different scenarios were used to estimate the hydraulic properties of C and T treatments (Table 3). They differ in their level of operationality for forest practitioners. The first scenario was based only on estimated values of hydraulic properties with the BEST method. It is named REF_BEST and is the least operational scenario as it requires complex measurements and algorithms to estimate hydraulic properties on both treatments. However, it is supposed to be the most precise and serves as reference. The two following scenarios are based on BEST-estimated values of the undisturbed soil at 0–10 cm depth (estimation of the $\mu + {R}_{s}$ term of Eqs. 3 and 4 for a new site) and on the values of hydraulic parameters estimated with SITE_mod (${\alpha }_{i}+{\beta }_{j}+{\alpha }_{i}:{\beta }_{j}$) and the PROX_mod (${\alpha }_{i}+{\gamma }_{l}+{\alpha }_{i}:{\gamma }_{l}$) models for the 15–25 cm layer of undisturbed soils (${\alpha }_{i}$) and both layers of the T-treatment soils (${\beta }_{j}+{\alpha }_{i}:{\beta }_{j}$ or ${\gamma }_{l}+{\alpha }_{i}:{\gamma }_{l}$). These two scenarios are named SITE and PROX scenarios and are also poorly operational. Indeed, they need the prior knowledge of soil hydraulic properties for the studied sites, which are rarely accessible to forest practitioners, which decreases the operational effectiveness of our predictive models. We then tested more operational scenarios that rely only on hydraulic properties obtained from textural properties and pedotransfer functions, PTF. We choose a class of PTF developed by Wessolek et al. (2009) from a German forest soils database to predict the hydraulic properties of the undisturbed soils. The advantages of this class of PTF are its validity for forest soils, its accessibility to forest managers, and its ease of use as it is based on soil textural class and provides the hydraulic parameters for the water retention and hydraulic conductivity curves. Then, the coefficients of SITE_mod and PROX_mod models were added to the PTF-estimated properties to estimate the properties of the skid trails: these two scenarios are named PTF_SITE and PTF_PROX respectively. These two scenarios are the most operational scenarios as they require textural class only.

Table 3 Description and summary of the data required to estimate soil hydraulic parameters of skid trails and level of operationality of the different approaches tested with water flow modelling

The hydraulic parameters estimated with these methods were used to model water flow in the different types of soils. The soil water flux was modelled with Richards’s equation solved using the finite difference approach. Soil profile was discretized into a grid composed of 80 nodes, distributed into three layers at 0–10, 10–30 cm and 30–100 cm. The hydraulic properties for the 0–10- and 10–30-cm layers were fixed according to the different scenarios (Table 3). The hydraulic properties for the 30–100-cm layer were equalled to those of the 10–30-cm layer, except for K_sat which was divided by 10. Simulations were performed over 3 years between 2016 and 2019 using SAFRAN daily climatic data. SAFRAN is a spatialization of observations from meteorological stations collected by Météo-France, the French national weather service.

All simulations with the different scenarios are based on the soil properties and forest stand data obtained on the AZZ25 site: roots distribution, Leaf Area Index (LAI), budburst and yellowing days. Simulated water contents were compared to measured water contents at three layers (10–20, 25–30 and 55–65 cm) for both treatments (Pousse et al. 2022). To compare simulated and observed water contents, daily means were computed for the two layers considered, i.e. 0–30 cm (mean of all the values between 0 and 30 cm) and 30–100 cm (mean of all the values between 31 and 100 cm). For each scenario, we computed the difference in simulated water saturation index between C and T treatments, hereafter named ∆compaction. We chose the saturation index instead of the water content to reduce bias induced by errors on the scaling parameters related to the water retention curves (residual and saturated water contents).

3 Results

3.1 Effect of compaction on soil hydraulic properties

The shapes of the retention curves at θ < θ_sat (i.e., the slope of the oblique part of the retention curves) were often the same between treatments (Fig. 4). Indeed, BEST methods use functions based only on particle size distributions for the estimation of the shape parameters (n, m and η) which determine the shapes of the water retention and the hydraulic conductivity curves. The variation of soil granulometry did not differ significantly between the C and T treatments over all the dataset (p-value = 0.84 for the clay content), leading to similar values of shape parameters. Therefore, the estimated values of the shape parameters and related shapes of retention curves were not expected to vary with soil compaction. The highest variation in clay content between treatments was found at ABB site. At 15–25 cm in plot ABB, we measured clay contents of 60 and 45% for respectively C and T treatments. The reason of this difference could be local heterogeneity and/or mixing soil layers during soil deformation. Soil volume loss induced by compaction or rutting can lead to horizons overlapping when comparing treatments, especially for soils with strong textural gradient (for example, a thin silty horizon on a clayey horizon).

Fig. 4

Mean retention curves per site representing the relation between the matrix potential h and the soil water content θ at two depths (continuous lines for 0–10 cm and dashed lines for 15–25 cm) and for two treatments (black line for the Control and red lines for the Trafficked treatment)

The scale parameter h_g is strongly related to the position of the inflexion points of the retention curves (null second derivative of h(θ)) and thus to the length of the plateaus close to saturation. High values of |h_g| indicate large plateaus and point at soils with very fine textures involving a majority of small pore sizes requiring large suctions to be desaturated (Lassabatere et al. 2021, 2023). For very fine soils, the beginning of desaturation may even require imposing a suction above a given value called the air-entry point and noted h_a. Such pattern is not included in the van Genuchten model implemented in BEST methods and should be taken into account when computing the sorptivity (Lassabatere et al. 2021, 2023). Treatment effect was not significant on h_g (i.e., p-value = 0.77, Fig. 5).

Fig. 5

Distribution of the measured parameters for the Control (C) and Trafficked (T) treatment at two depths. K_sat, θ_sat and h_g are respectively the hydraulic conductivity, the water content at saturation and the value of the matric potential at the inflexion point of the retention curve

Variations of the soil water content at saturation, θ_sat, define the different positions of the vertical plateaus close to saturation (Fig. 4). The interaction between treatment and depth was significant on θ_sat. On average, the values of θ_sat are reduced by 0.05 m³ m⁻³ on trafficked soils at 15–25 cm and by 0.02 m³ m⁻³ at 0–10 cm depth (Fig. 5, Table 4). Logically, soil compaction significantly reduces soil porosity and saturated water content. For some sites, the difference in θ_sat between C- and T-treatments at 0–10 cm was small to negligible, i.e. sites ABB, ARF2, AZ25, CAM, H, and VER11 (Fig. 6). The values on the trafficked treatment of the two monitored experimental sites (AZ25 and CA, “experimental” in Fig. 7) did not significantly differ from the values of the management sites where infiltration tests were realized 1 to 2 years before the upcoming harvest without monitoring of the traffic intensity (“management” in Fig. 7). The absence of outliers was detected with the analysis of the standardized marginal residuals of the SITE_mod models. These two monitored experimental sites behaved like the “forest management” sites.

Table 4 Relationship between treatment, depth and soil hydraulic parameters while considering the site as random effect (= SITE_mod)

Fig. 6

Distribution of the measured soil water content at saturation at 0–10 cm as a function of treatment and site

Fig. 7

Measured versus predicted values of the T-treatment with SITE_mod model for θ_sat and log₁₀(K_sat) in A aborted experiment exclusion case and B aborted experiment inclusion case with K_sat fixed at 0.0001 mm s⁻¹ for aborted experiment. Grey points are the measured and predicted values for the two experimental sites where traffic intensity and time since traffic is controlled (AZ25 and CA plots), black points are the measured and predicted values for the remaining plots where traffic records and details on machinery used is unknown

Traffic significantly reduced the hydraulic conductivity at saturation (K_sat), by dividing the conductivity of the C-treatment by 2.4 for all soil layers (Fig. 5, Table 4).

We recall that, for some sites, the saturated hydraulic conductivities could not be determined with the BEST methods because of too small cumulative infiltrations. In the following, we test the abovementioned strategies to replace these unknown values. When incorporating aborted experiments in the dataset considering scenarios 2–4 (fixed K_sat values assigned to aborted experiments), we increased the variability of K_sat leading to RMSE of log₁₀(K_sat) between 0.50 and 0.75 when K_sat of aborted tests was changed from K_sat = 10⁻⁴ to K_sat = 10⁻⁶ mm s⁻¹ (Table 4). Then, the coefficient of determination decreased from 0.55 to 0.47. Hereafter, we focused on the data set in which aborted experiment were removed (hereafter referred to as the aborted experiment exclusion case) and on the data set including aborted experiment with fixed K_sat value leading the lowest model error and the highest R², i.e. K_sat = 10⁻⁴ mm s⁻¹ (hereafter referred to as the aborted experiment inclusion case). The interaction between treatment and site effects becomes significant in the aborted experiment inclusion case and reflects the weight of these trafficked low-permeability soils in the assessment of the T-treatment effect. The K_sat values of these low-permeability soils were mostly poorly predicted (Fig. 7B, measured log₁₀(K_sat) = − 4), probably because Beerkan experiments abortion prevented to estimate very low values of K_sat with the BEST methods (Fig. 7A, minimum measured log₁₀(K_sat) > − 4).

3.2 Effect of the intensity of skid trail degradation

Several types of degradation scores were tested to predict the observed impact on θ_sat and log₁₀(K_sat): rut depth, rut depth class (≤ 5 cm: 16 sites, 5–10 cm: 19 sites, > 10 cm: 6 sites), a qualitative score used by Frey et al. (2009) including rut depth class and change in soil structure (score 1: 3 sites, score 2: 29 sites, score3: 9 sites), and our score proposed in this study including rut depth class, change in soil structure and depth of redoximorphic features (Fig. 2, score 1: 6 sites, score 2: 23 sites, score 3: 12 sites). The description of the intensity of skid trail degradation with our score gave the best results to predict the observed impact on θ_sat and log₁₀(K_sat). Hereafter, we focused on the implementation of this qualitative score in the PROX_mod model (Table 5, Fig. 8). In addition, the focus is put on the analysis of θ_sat and log₁₀(K_sat) since the other variables proved to be less sensitive to soil compaction.

Table 5 Relationship between the intensity of skid trail degradation, depth and soil hydraulic parameters while considering the site as random effect (= PROX_mod)

Fig. 8

Measured versus predicted values of the T-treatment with PROX_mod model for θ_sat and log₁₀(K_sat) in A aborted experiment exclusion case and B aborted experiment inclusion case with K_sat fixed at 0.0001 mm s⁻¹ for aborted experiment. Grey points are the measured and predicted values for the two experimental sites where traffic intensity and time since traffic is controlled (AZ25 and CA plots), black points are the measured and predicted values for the remaining plots where traffic records and details on machinery used is unknown

For PROX_mod, both depth and intensity of skid degradation had an impact on θ_sat and log₁₀(K_sat). In addition, their interactions depth*intensity of skid trail degradation were significant for both θ_sat and log₁₀(K_sat). The effect on θ_sat at 0–10-cm increases from soil degradation score 1 to 2 and decreases for the soil degradation score 3 (Table 5). At 15–25 cm, θ_sat decreases when the soil degradation score increases. For log₁₀(K_sat), the effect decreases when the soil degradation score increases at 0–10 cm. At 15–25 cm, the effect increases from score degradation 1 to 2 but decreases from 2 to 3. The effect is higher when the aborted experiments are included (aborted experiments inclusion versus aborted experiments exclusion cases). PROX_mod gave better predictions than SITE_mod, only for K_sat (Tables 4 and 5, Figs. 7 and 8).

3.3 Use of predicted soil hydraulic properties in a water flow model

Predicted θ_sat for the T-treatment, using the PTF of Wessolek et al. (2009) and SITE_mod or PROX_mod coefficients, are always less than the measured values, except for one point, and predicted K_sat are mainly less than the reference values of K_sat estimated with BEST methods (i.e. example for the PTF_SITE scenario in Fig. 9). We then expected errors in terms of modelling water flow with the use of those hydraulic parameters.

Fig. 9

Measured versus predicted values of the T-treatment for the PTF_SITE scenario; using the PTF developed by Wessolek et al. (2009) to predict the hydraulic properties of the undisturbed soils and the SITE_mod model to predict θ_sat and log₁₀(K_sat) of the T-treatment in A aborted experiment exclusion case and B aborted experiment inclusion case with K_sat fixed at 0.0001 mm s⁻¹ for aborted experiment

The difference in water saturation index between C and T treatments, ∆compaction, obtained with scenarios SITE, PROX, PTF_SITE and PTF_PROX (Table 3) are compared to ∆compaction calculated with the REF_BEST scenario and to the observed data (Fig. 10). For the REF_BEST scenario, hydraulic properties were estimated from infiltration experiment with BEST method for both C and T treatments. For the SITE or PROX scenarios, hydraulic properties were estimated from BEST-estimated properties for the C-treatment and from their incorporation into the SITE_mod or PROX_mod models for the T-treatment. For the PTF_SITE and PTF_PROX scenarios, hydraulic properties were estimated from the PTF developed by Wessolek et al. (2009) for the C-treatment and then incorporated into SITE_mod or PROX_mod models to derive those of the T-treatment.

Fig. 10

Difference in simulated water saturation index between C and T treatments (y-axis) by a water flow model parametrized with hydraulic properties estimated from the different scenarios described in Table 3 for the AZ25 site. Saturation index of the C-treatment is higher than the one of the T-treatment when the curves are above the black dotted line. OBS: observed differences in saturation index between treatments, REF_BEST: hydraulic properties estimated from infiltration experiment with BEST method for both treatments, SITE: hydraulic properties estimated from infiltration experiment with the BEST method for the C-treatment and from the values of the C-treatment and SITE_mod model for the T-treatment, PROX: hydraulic properties estimated from infiltration experiment with the BEST method for the C-treatment and from the values of the C-treatment and PROX_mod model for the T-treatment, PTF_SITE: hydraulic properties for the T-treatment estimated from the PTF developed by Wessolek et al. (2009) for the C-treatment and SITE_mod model, PROX: hydraulic properties for the T-treatment estimated the PTF developed by Wessolek et al. (2009) for the C-treatment and PROX_mod model

∆compaction from REF_BEST scenario was always negative for all soil layers meaning that soil saturation index of the T-treatment was always higher than the one of the C-treatment. The measured soil saturation index of the T-treatment was also always higher than the one of the C-treatment for the 10 to 30-cm layer but not for the 30 to 100-cm layer (Fig. 10 top versus bottom). Besides the ∆compaction from REF_BEST scenario was higher than the observed one. The root mean square error, RMSE, between observed ∆compaction and the REF_BEST one is high compared to other scenarios (Table 6). Temporal variations of the ∆compaction for the REF_BEST scenario and the observations were similar, except for some discrepancies that may originate from divergences in climatic data between SAFRAN and local weather and in hydraulic properties as BEST method underestimate treatment effect on them.

Table 6 RMSE between observed and predicted differences in simulated water saturation index between C and T treatments (Fig. 10), according to the approaches described in Table 3

Evolutions over time of ∆compaction simulated by the four other scenarios are similar to the ones obtained with the observations. The order of magnitude of ∆compaction is also comparable except for the PTF_PROX scenario. Results of the PTF_PROX scenario were the farthest to observations (Fig. 10) even if the RMSE over the entire period was the same as for the other scenarios (Table 6). ∆compaction evolved only slightly as a function of time and depth and its magnitude was very low compared to the other scenarios. The PTF_SITE, SITE and PROX scenarios were close to each other and gave better results for the 10–30-cm layer than for the 30–100-cm layer. We expected that the results obtained with SITE and PROX scenarios would be closer to the results of the REF_BEST scenario since hydraulic properties were shared for the C treatment. It is not always the case with periods where the data from the PTF_SITE are closer to the REF_BEST. The PTF_SITE scenario is operational and, compared to the SITE, PROX or REF_BEST scenarios, does not decrease accuracy for the implementation of decision support tools based on the use of water flow models, without requiring precise knowledge of past harvesting conditions nor time-consuming measurements of cumulative infiltrations.

4 Discussion

4.1 Estimation of soil hydraulic parameters and effect of compaction

Our results showed that a compaction effect was only detected for θ_sat and log₁₀(K_sat). We expected that the scale parameter h_g of the water retention curve would be modified for compacted soils since it is related to the average pore size which is diminished by soil compaction. In addition, it may be considered as a proxy of the air-entry point of the water retention curve that is expected to increase in absolute value with a decrease in macroporosity (Dexter 2004). However, the compaction effect on this parameter was not significant in our study. Working on agricultural soils, Richard et al. (2001) showed that the pore size distribution is affected by compaction especially for the largest pores: after compaction of a silty agricultural soil, they observed a loss of 20% of total porosity and a loss of 70% of pores of size between 40 and 360 µm. Bottinelli et al. (2014) showed that the specific surface porosity of large pores was considerably reduced after compaction for the AZ25 and CA experimental sites. On these two experimental sites, pore size distribution was impacted by traffic (Goutal-Pousse et al. 2016), whereas no impact on h_g was evidenced by our study. As h_g is the last parameter estimated by the BEST methods, its estimation depends on the values of all variables measured in the field (θ_initial, θ_sat) and all parameters estimated previously by the BEST method (n, η, S, K_sat). Consequently, we might expect a larger uncertainty for the estimation of this parameter than for the others. This may explain why there are no clear trends between treatments and this parameter seems to be insensitive to the compaction effect.

Shape parameters of the water retention curves are expected to vary as compaction modifies the pore size distribution by reducing macroporosity and by modifying the pore shape. Some authors (e.g. Dexter 2004; Naderi-Boldaji and Keller 2016) used the slope of the retention curve at the inflexion point as a soil physical quality indicator for compaction since its value decreases with the increase in bulk density according to their findings. The derivation of the slope at the inflexion point, based on the van genuchten model for the water retention curve, reveals that this slope is a function of the parameter n (Dexter 2004). The use of the bulk density in addition to particle size distribution in the BEST functions might improve the assessment of compaction effects on n and, consequently, h_g. This gives room to further research on BEST methods and the ways of combining textural and structural analysis in these methods.

Aborted tests, where infiltration time was too long to be measured on the field, only occurred in the T-treatment and were ascribed to very low values of soil hydraulic conductivity. The proposed method to fix values for K_sat for aborted tests slightly improved the assessment of the effect of compaction on K_sat. The calculated compaction effect was higher, but the modelling error also increased. Besides, BEST simulation errors were high in the T-treatment compared to the C-treatment (Martin et al. 2024) and two sites displaying severe degradation classes (BTG and FTG) had to be excluded from the dataset because of high BEST simulation errors. These high BEST simulation errors were related to the high concavity of the cumulated infiltration curves, particularly at the beginning of the curves (Martin et al. 2024). The high concavity of the infiltration curves in the T-treatment may be associated with soil sealing (Di Prima et al. 2018) that could be due to horizontal pores which were observed in the T-treatment at the AZ25 and CA sites (Bottinelli et al. 2014). Therefore, the calculated effect of traffic and related compaction on K_sat might be underestimated, especially for severe degradation classes.

4.2 Role of soil degradation intensity and soil evolution since the last degradation

With BEST methods, total porosity and θ_sat are calculated in the same way using soil bulk density and solid (or mineral) density. The low impact of traffic on θ_sat at 0–10 cm observed in several sites could be due to soil restoration between two forest operations, because some of them still had high visual degradation scores (including rut depth, change in structure, and apparition of redoximorphic features) and a high impact on θ_sat at 15–25 cm. Besides, the AZ25 and CA sites, where initial impact and soil evolution since initial impact are monitored, had the same behaviour than the rest of the sites. In both sites, the high initial impact on total porosity between 0 and 7 to 10 cm decreased with time, whereas the initial impact observed below 7–10 cm did not evolve and remained significant (Goutal et al. 2012; Bottinelli et al. 2014). Other studies on forest soils pointed at the same short-term (1 to 10 years) evolution in the first 10 cm (Fründ and Averdiek 2016) and a faster evolution on clayey and nutrients’ rich soils than on sandy and poor soils (Ebeling et al. 2017; Mohieddinne et al. 2019). We hypothesised that soil pH is a relevant parameter for explaining differences in soil recovery dynamics, especially the threshold pH of 4.5 which is also known to be important for aluminium speciation and soil biogeochemical functioning (Driscoll and Schecher 1990). Besides, a soil pH threshold of 4.5 explained the varying evolution of soil solution after traffic at AZ25 and CA sites (Ranger et al. 2021). In our study, soil pH was linked to the difference in θ_sat between the C and T treatments at 0–10 cm; i.e. the difference was high for soils displaying pH under 4.5 and low for soils displaying pH above 4.5, except for the CAM site which is the sandiest site of our dataset (Fig. 11). However, more data are needed in order to confirm this relationship between soil acidity and soil porosity restoration dynamics at 0–10 cm.

Fig. 11

Percentage of reduction between C and T treatments for measured θ_sat as a function of soil pH in water at 0–10 cm

The use a visual score of skid trail degradation intensity did not improve the estimation of θ_sat for skid trails, yet it gave trends worth further investigations. Besides it did improve the estimation of K_sat. Frey et al. (2009) have shown that the strongest effects of compaction concerned severe degradation classes whereas Aust et al. (1998) have shown that the effect of compaction on bulk density and saturated hydraulic conductivity is the same between rutted and deeply rutted forest soils. Aust et al. (1998) have explained this fact by the incompressibility of water when soil porosity was saturated at the time of traffic. This may explain the decrease of the effect of disturbance on θ_sat at 0–10 cm and on K_sat between score 2 and 3 (PROX_mod).

For degradation score 1, the impact on saturated hydraulic conductivity was high at 0–10 cm and decreased with depth. For degradation scores 2 and 3, the impact on θ_sat and K_sat were higher at 15–25 cm than at 0–10 cm, especially for degradation score 3 where the increase in impact from 0–10 cm to 15–25 cm was the highest. This points out that the measured impacts for score degradation 3 could increase at depths deeper than 25 cm as soil stresses reach deeper depths with increasing traffic degradation intensity, i.e. increasing water content at the time of traffic (Smith et al. 2000) and/or increasing machinery weight (Keller et al. 2019). Yet, measurements done at 30–40 cm depth were mostly unsuccessful and we were unable to perform the analysis of soil compaction effect for deep soil layers using the Beerkan protocol and BEST methods.

4.3 Comparison of the different approaches used to predict soil hydraulic parameters

The use of mechanistic model of water flow is an approach to support decision-making process. In our study the hydraulic parameters h_g, n and η were derived from the undisturbed soils (estimated with the BEST method or a PTF) while K_sat and θ_sat were estimated using an estimated compaction effect (SITE_mod or PROX_mod). The PTF_SITE scenario used in a mechanistic model of water flow could be an operational way to predict soil saturation index of the first 30 cm and help forest managers to better plan skidding operations. Yet, the proposed method to infer soil hydraulic properties using the outcome of our study should be tested on more sites and with several PTFs. Besides, the construction of hydraulic functions with parameters for the compacted soil coming from different origins might not be optimal (Wessolek et al. 2009).

5 Conclusions

The proposed database and related statistical analysis did not include very fine textured soils and is limited to rut having a depth below 12 cm, i.e. , it does not include severe traffic-induced degradation. The validity of our conclusions is therefore limited to medium to coarse textured soils with low level of traffic-induced degradation. We expect a larger range of variations in θ_sat, K_sat and the other hydraulic parameters when enlarging the range of soil and compaction conditions. Therefore, this first attempt at predicting the hydraulic properties of compacted forest soils must be completed by other measurements to oversee the entire range of textures and levels of degradation. Adding more measurements on skid trails with ruts above 12 cm would reinforce our database to validate or invalidate the detected trends between hydraulic properties and visual score of skid trail degradation, i.e., increasing depth of maximal impact with increasing degradation score and low impact on bulk density associated with high impact on K_sat for the highest degradation score. However, treating more severely degraded skid trails with deeper ruts seems difficult to realize with the BEST procedure, because of the high risk of failure of infiltration measurement or of poor fit of the BEST algorithms. Alternative methods should be searched. Moreover, the effects of compaction on the scale parameter of the water retention curve should be investigated.

Our approach has validated the following strategy for the hydraulic characterization of forest plots with skid trails. The hydraulic parameters of undisturbed soils can be measured with infiltration experiments or estimated from a PTF that is the most operational option for forest practitioners. Then, the hydraulic properties of compacted soils under skid trails may be predicted accurately with our proposed SITE_mod model. In addition, the use of those predicted properties proved efficient for the modelling of water flow in the soil profile. Indeed, the saturation degree with numerical modelling based on the estimated hydraulic parameters proved close to the observed values for a specific soil profile. However, some discrepancies were observed during some periods of the year and in deeper soil layers and further work is still needed to reconstruct the hydraulic functions for the T-treatment (compacted and rutted soils) based on undisturbed soils. Nevertheless, our simulations highlight the change in soil water content dynamics between compacted and uncompacted soils. This is encouraging for the implementation of future decision support tools based on the use of water flow models, without requiring precise knowledge of past harvesting conditions or soil degradation status.