<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

Manon Martin¹ ID, André Chanzy¹ ID, Laurent Lassabatere² ID, Arnaud Legout³ ID, Noémie Pousse⁴ ID and Stéphane Ruy¹* ID

**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_{10}(K_{sat})$) and the soil water content at saturation ($\theta_{sat}$). For the skid trails, $\theta_{sat}$ was reduced by - 0.02 and -0.11 $m^3 m^{-3}$ in the 0–10 cm and 15–25 cm layers respectively, compared to undisturbed topsoil (0–10 cm). $\log_{10}(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.

<figure><img src="image_4.png" /></figure>
**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).

**Keywords** Soil hydraulic parameters, Compaction, Forest soil, Pedotransfer function, BEST method

# 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

<table><caption>**Table 1** Name, location, stand composition of the studied plots and intensity of degradation of the skid trails</caption><thead><tr><th>Plot name</th><th>GPS coordinates (decimal degrees) lat, long</th><th>Alt (m)</th><th>Stand</th><th>Age (year)</th><th>Degradation class of the skid trail</th><th>Rut depth (cm)</th></tr></thead><tbody><tr><td>ABB</td><td>46.80, 2.60</td><td>170</td><td>Oak</td><td>70</td><td>Moderate</td><td>7.5</td></tr><tr><td>ARF1</td><td>43.40, 2.20</td><td>700</td><td>Beech</td><td>45</td><td>Moderate</td><td>7</td></tr><tr><td>ARF2</td><td>43.40, 2.20</td><td>700</td><td>Beech</td><td>45</td><td>Moderate</td><td>12</td></tr><tr><td>AZ</td><td>48.51, 6.69</td><td>326</td><td>Beech</td><td>90</td><td>Low</td><td>3</td></tr><tr><td>AZ25</td><td>48.50, 6.70</td><td>320</td><td>Oak</td><td>15</td><td>Severe</td><td>5</td></tr><tr><td>BTG</td><td>48.82, 6.98</td><td>256</td><td>Oak, beech</td><td>35</td><td>Severe</td><td>7.5</td></tr><tr><td>CA</td><td>49.14, 5.02</td><td>272</td><td>Oak</td><td>15</td><td>Severe</td><td>5</td></tr><tr><td>CAM</td><td>43.64, 2.92</td><td>1000</td><td>Spruce</td><td>60</td><td>Low</td><td>4</td></tr><tr><td>FTG</td><td>48.82, 6.91</td><td>240</td><td>Oak, hornbeam, ash</td><td>50</td><td>Severe</td><td>15</td></tr><tr><td>FTG148</td><td>48.80, 6.96</td><td>246</td><td>Oak, hornbeam</td><td>55</td><td>Low</td><td>3</td></tr><tr><td>H</td><td>48.67, 7.06</td><td>306</td><td>Beech, hornbeam</td><td>60</td><td>Severe</td><td>10</td></tr><tr><td>POC1</td><td>47.08, 5.32</td><td>186</td><td>Oak</td><td>70</td><td>Moderate</td><td>1</td></tr><tr><td>POC2</td><td>47.08, 5.32</td><td>186</td><td>Oak</td><td>70</td><td>Moderate</td><td>7</td></tr><tr><td>POU2</td><td>46.98, 5.24</td><td>194</td><td>Oak, beech</td><td>80–100</td><td>Moderate</td><td>4.5</td></tr><tr><td>POU3</td><td>46.98, 5.24</td><td>194</td><td>Oak, beech</td><td>80–100</td><td>Moderate</td><td>7.5</td></tr><tr><td>SAU</td><td>47.47, 5.73</td><td>244</td><td>Oak, birch</td><td>25</td><td>Severe</td><td>7</td></tr><tr><td>SRG</td><td>48.77, 6.95</td><td>268</td><td>Spruce, oak</td><td>40</td><td>Severe</td><td>17</td></tr><tr><td>VER11</td><td>47.20, 6.28</td><td>578</td><td>Beech, oak</td><td>50–60</td><td>Moderate</td><td>9</td></tr><tr><td>VER6</td><td>47.20, 6.28</td><td>592</td><td>Beech</td><td>50–60</td><td>Moderate</td><td>12</td></tr></tbody></table>

Ash: *Fraxinus excelsior* L., beech: *Fagus sylvatica* L., birch *Betula pendula* Roth, hornbeam: *Carpinus betulus* L., oak: *Quercus petraea* Liebl. and *Quercus robur* L., spruce: *Picea abies* (L.) H. Karst

Degradation class of the skid trail is determined according to rut depth, soil structure and the apparition depth of redoximorphic features (see Fig. 2): 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. 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
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Location of studied plots</figcaption></figure>

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

## 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):

* 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 ($\rho_b$) and the initial water content $\theta_0$. 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).

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

<figure><img src="image_6.png" /><figcaption>**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</figcaption></figure>
<table><caption>**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)</caption><thead><tr><th>Site name</th><th>Depth (cm)</th><th>N</th><th>Clay (%)</th><th>Silt (%)</th><th>C<sub>org</sub> (g kg<sup>-1</sup>)</th><th>pH<sub>water</sub></th><th>USDA texture</th><th>HYPRES texture</th></tr></thead><tbody><tr><td rowspan="2">ABB</td><td>0-10</td><td>2</td><td>41.3</td><td>41.6</td><td>58.9</td><td>5.4</td><td>SiCl</td><td>F</td></tr><tr><td>15-25</td><td>2</td><td>52.4</td><td>34.0</td><td>22.8</td><td>5.3</td><td>Cl</td><td>F</td></tr><tr><td rowspan="2">ARF1</td><td>0-10</td><td>2</td><td>20.8</td><td>23.0</td><td>72.6</td><td>4.4</td><td>SaClLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>16.7</td><td>28.6</td><td>46.7</td><td>4.4</td><td>SaLo</td><td>M</td></tr><tr><td rowspan="2">ARF2</td><td>0-10</td><td>2</td><td>20.8</td><td>26.9</td><td>59.0</td><td>5.0</td><td>SaClLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>18.0</td><td>31.8</td><td>40.0</td><td>4.7</td><td>Lo</td><td>M</td></tr><tr><td rowspan="2">AZ</td><td>0-10</td><td>2</td><td>14.5</td><td>66.5</td><td>42.0</td><td>4.2</td><td>SiLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>15.3</td><td>66.0</td><td>8.9</td><td>4.2</td><td>SiLo</td><td>M</td></tr><tr><td rowspan="2">AZ25</td><td>0-10</td><td>30</td><td>19.9</td><td>57.9</td><td>25.2</td><td>4.8</td><td>SiLo</td><td>M</td></tr><tr><td>15-25</td><td>30</td><td>19.4</td><td>58.8</td><td>10.3</td><td>4.6</td><td>SiLo</td><td>M</td></tr><tr><td>BTG</td><td>0-10</td><td>2</td><td>18.2</td><td>70.2</td><td>-</td><td>-</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">CA</td><td>0-10</td><td>30</td><td>11.9</td><td>74.2</td><td>28.3</td><td>4.4</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>30</td><td>13.4</td><td>72.4</td><td>7.2</td><td>4.5</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">CAM</td><td>0-10</td><td>2</td><td>18.8</td><td>16.2</td><td>115.0</td><td>3.9</td><td>SaLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>11.3</td><td>16.2</td><td>45.9</td><td>4.3</td><td>SaLo</td><td>C</td></tr><tr><td>FTG</td><td>0-10</td><td>2</td><td>26.6</td><td>59.6</td><td>-</td><td>-</td><td>SiClLo</td><td>MF</td></tr><tr><td>FTG148</td><td>0-10</td><td>2</td><td>16.4</td><td>63.7</td><td>-</td><td>-</td><td>SiLo</td><td>M</td></tr><tr><td rowspan="2">H</td><td>0-10</td><td>2</td><td>22.9</td><td>62.4</td><td>60.3</td><td>6.1</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>2</td><td>24.1</td><td>61.8</td><td>17.6</td><td>6.9</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">POC1</td><td>0-10</td><td>2</td><td>14.2</td><td>76.9</td><td>39.3</td><td>4.4</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>2</td><td>14.5</td><td>76.9</td><td>15.7</td><td>4.4</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">POC2</td><td>0-10</td><td>2</td><td>12.8</td><td>67.8</td><td>10.9</td><td>4.5</td><td>SiLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>16.4</td><td>64.7</td><td>4.5</td><td>4.7</td><td>SiLo</td><td>M</td></tr><tr><td rowspan="2">POU2</td><td>0-10</td><td>2</td><td>17.2</td><td>70.7</td><td>144.3</td><td>3.8</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>2</td><td>15.0</td><td>72.8</td><td>24.3</td><td>4.1</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">POU3</td><td>0-10</td><td>2</td><td>17.6</td><td>55.7</td><td>24.4</td><td>4.6</td><td>SiLo</td><td>M</td></tr><tr><td>15-25</td><td>2</td><td>20.0</td><td>54.0</td><td>13.2</td><td>4.5</td><td>SiLo</td><td>M</td></tr><tr><td>SAU</td><td>0-25</td><td>1</td><td>12.7</td><td>52.8</td><td>-</td><td>-</td><td>SiLo</td><td>M</td></tr><tr><td>SRG</td><td>0-10</td><td>2</td><td>37.0</td><td>55.0</td><td>-</td><td>-</td><td>SiClLo</td><td>F</td></tr><tr><td rowspan="2">VER11</td><td>0-10</td><td>2</td><td>22.1</td><td>63.2</td><td>16.9</td><td>4.8</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>2</td><td>26.0</td><td>61.8</td><td>32.6</td><td>4.5</td><td>SiLo</td><td>MF</td></tr><tr><td rowspan="2">VER6</td><td>0-10</td><td>1</td><td>27.8</td><td>65.0</td><td>-</td><td>4.6</td><td>SiLo</td><td>MF</td></tr><tr><td>15-25</td><td>1</td><td>27.4</td><td>66.9</td><td>5.9</td><td>4.7</td><td>SiLo</td><td>MF</td></tr></tbody></table>

Cl Clay, SiCl Silty clay, SiLo Silt loam, SiClLo Silty clay loam, Lo Loam, SaLo Sandy loam, SaClLo Sandy clay loam, VF Very Fine, F Fine, MF Medium Fine, M Medium, C Coarse. N number of samples

(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} \quad (1)$$

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

where $\theta_{sat}$ and $\theta_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). $\theta_{sat}$ is deduced from the bulk density $\rho_b$ measurement and the particle density ($\rho_s$=2.65 g cm⁻³): $\theta_{sat}=1 - \rho_b / \rho_s$ and $\theta_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^{-4} \ mm \ s^{-1}$, all experiments should have led to success whereas when $K_{sat} < 10^{-4} \ mm \ s^{-1}$ all experiments should have failed. This range was therefore considered relevant for the definition of the majorant of $K_{sat}$ for aborted experiments.

<figure><img src="image_7.png" /><figcaption>**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</figcaption></figure>
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 Ksat for all aborted experiments. Values of Ksat 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ₗ 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₃ 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₃ₗ 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.

$$ \hat{y}_{ijk} = \mu + \alpha_i + \beta_j + \alpha_i : \beta_j + R_s + R_{sj} + \varepsilon_{ijk} \quad (3) $$

where $\hat{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, αᵢ and βⱼ are the fixed effects with αᵢ the depth effect (αᵢ= 0 at 0–10 cm) and βⱼ the compaction effect (βⱼ= 0 for C-treatment), αᵢ : βⱼ is the interaction effect (traffic effect at 15–25 cm depth if significantly different from μ + αᵢ + βⱼ), R₃ is the site random effect, R₃ₗ is the site random effect for treatment j and ε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.

$$ \hat{y}_{ilk} = \mu + \alpha_i + \gamma_l + \alpha_i : \gamma_l + R_s + \varepsilon_{ilk} \quad (4) $$

where $\hat{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, αᵢ is the depth effect (αᵢ= 0 at 0–10 cm), γᵢ is the coefficient for the level l of visual indicator of skid trail degradation (γᵢ=0 for the C-treatment), αᵢ : γᵢ is the interaction effect between depth and visual indicator of skid trail degradation, R₃ is the site random effect (related to the between-site variability of soil properties for the C-treatment at 0–10 cm) and ε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_practicalability_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 hydrau-
lic properties of C and T treatments (Table 3). They
differ in their level of operationality for forest practi-
tioners. 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 sce-
narios are based on BEST-estimated values of the undis-
turbed soil at 0–10 cm depth (estimation of the µ + R_s
term of Eqs. 3 and 4 for a new site) and on the values
of hydraulic parameters estimated with SITE_mod
(α_i + β_j + α_i : β_j) and the PROX_mod (α_i + γ_l + α_i : γ_l)
models for the 15–25 cm layer of undisturbed soils (α_i)
and both layers of the T-treatment soils (β_j + α_i : β_j
or γ_l + α_i : γ_l). These two scenarios are named SITE
and PROX scenarios and are also poorly operational.
Indeed, they need the prior knowledge of soil hydrau-
lic 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 prop-
erties 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 hydrau-
lic 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 hydrau-
lic 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.

The hydraulic parameters estimated with these meth-
ods 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 differ-
ent scenarios (Table 3). The hydraulic properties for the
30–100-cm layer were equalled to those of the 10–30-cm

layer, except for Ksat which was divided by 10. Simula-
tions were performed over 3 years between 2016 and
2019 using SAFRAN daily climatic data. SAFRAN is a
spatialization of observations from meteorological sta-
tions 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 lay-
ers (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 distri-
butions 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 high-
est variation in clay content between treatments was
found at ABB site. At 15–25 cm in plot ABB, we meas-
ured 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 defor-
mation. Soil volume loss induced by compaction or rut-
ting can lead to horizons overlapping when comparing
treatments, especially for soils with strong textural gradi-
ent (for example, a thin silty horizon on a clayey horizon).

The scale parameter *h*<sub>*g*</sub> is strongly related to the posi-
tion 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*<sub>*g*</sub>| indicate
large plateaus and point at soils with very fine textures
<table><caption>**Table 3** Description and summary of the data required to estimate soil hydraulic parameters of skid trails and level of operatibility of the different approaches tested with water flow modelling</caption><thead><tr><th>Approaches</th><th>Equation</th><th>Hydraulic parameters on Control soils</th><th>Hydraulic parameters on Trafficked soils</th><th>Summary of the needed data</th><th>Operatibility</th></tr><tr><th></th><th></th><th></th><th>Significant to compaction</th><th></th><th></th></tr></thead><tbody><tr><td>REF_BEST</td><td></td><td>K<sub>sat</sub>, θ<sub>sat</sub>, n, η and h<sub>g</sub> parameters are estimated from infiltration experiment on Control soils</td><td>K<sub>sat</sub>, θ<sub>sat</sub>, n, η and h<sub>g</sub> are estimated from infiltration experiment on Trafficked soils</td><td>Infiltration experiment with BEST method on Control and Trafficked soils</td><td>Low</td></tr><tr><td>SITE</td><td>y<sub>kij</sub> = μ<sub>measured</sub> + α<sub>i</sub> + β<sub>j</sub></td><td></td><td>Coefficients of SITE_mod are applied to log<sub>10</sub>(K<sub>sat</sub>) and θ<sub>sat</sub> of Control soils at 0–10 cm (depending on treatment and depth)</td><td>Infiltration experiment with BEST method on Control topsoil</td><td>Moderate</td></tr><tr><td>PROX</td><td>y<sub>kil</sub> = μ<sub>measured</sub> + α<sub>i</sub> + γ</td><td></td><td>Coefficients of PROX_mod are applied on log<sub>10</sub>(K<sub>sat</sub>) and θ<sub>sat</sub> of Control soils at 0–10 cm (depending on soil disturbance score and depth)</td><td>Infiltration experiment with BEST method on Control topsoil Soil disturbance classification</td><td>Moderate</td></tr><tr><td>PTF_SITE</td><td>y<sub>kij</sub> = μ<sub>PTF</sub> + α<sub>i</sub> + β<sub>j</sub></td><td>K<sub>sat</sub>, θ<sub>sat</sub>, n, η and h<sub>g</sub> parameters are estimated using a PTF applied to Control soils</td><td>Coefficients of SITE_mod are applied on log<sub>10</sub>(K<sub>sat</sub>) and θ<sub>sat</sub> of Control soils at 0–10 cm (depending on treatment and depth)</td><td>Soil textural class of Control soils (as PTF input)</td><td>High</td></tr><tr><td>PTF_PROX</td><td>y<sub>kol</sub> = μ<sub>PTF</sub> + α<sub>i</sub> + γ</td><td></td><td>Coefficients of PROX_mod are applied on log<sub>10</sub>(K<sub>sat</sub>) and θ<sub>sat</sub> of Control soils at 0–10 cm (depending on soil disturbance score and depth)</td><td>Soil textural class of Control soils (as PTF input) Soil disturbance classification</td><td>High</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**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)</figcaption></figure>

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*<sub>a</sub>. 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*<sub>g</sub> (i.e.,
*p*-value = 0.77, Fig. 5).
Variations of the soil water content at saturation, $\theta_{sat}$, define the different positions of the vertical plateaus close to saturation (Fig. 4). The interaction between treatment and depth was significant on $\theta_{sat}$. On average, the values of $\theta_{sat}$ are reduced by 0.05 m$^3$ m$^{-3}$ on trafficked soils at 15–25 cm and by 0.02 m$^3$ m$^{-3}$ 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 $\theta_{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

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Distribution of the measured parameters for the Control (C) and Trafficked (T) treatment at two depths. $K_{sat}$, $\theta_{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</figcaption></figure>

<table><caption>**Table 4** Relationship between treatment, depth and soil hydraulic parameters while considering the site as random effect (= SITE\_mod)</caption><thead><tr><th></th><th>θ<sub>sat</sub><sup>A</sup><br>m<sup>3</sup> m<sup>-3</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>A</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>B</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>C</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>D</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th></tr></thead><tbody><tr><td>μ (C at 0–10 cm)</td><td>0.64 (±0.006)</td><td>-1.72 (±0.05)</td><td>-1.74 (±0.05)</td><td>-1.74 (±0.05)</td><td>-1.75 (±0.06)</td></tr><tr><td>β<sub>j</sub> (T)</td><td>-0.019 (±0.010)</td><td>-0.38 (±0.05)</td><td>-0.50 (±0.06)</td><td>-0.60 (±0.06)</td><td>-0.69 (±0.08)</td></tr><tr><td>α<sub>1</sub> (15–25 cm)</td><td>-0.058 (±0.009)</td><td>-0.47 (±0.05)</td><td>-0.47 (±0.06)</td><td>-0.47 (±0.06)</td><td>-0.50 (±0.09)</td></tr><tr><td>α<sub>2</sub>β<sub>j</sub> (T at 15–25 cm)</td><td>-0.035 (±0.015)</td><td>p-value = 0.51</td><td>p-value = 0.57</td><td>p-value = 0.50</td><td>p-value = 0.47</td></tr><tr><td>Variance of R<sub>s</sub> (Site effect)</td><td>0.002</td><td>0.13</td><td>0.16</td><td>0.14</td><td>0.011</td></tr><tr><td>p-value of R<sub>sj</sub> (interaction between Treatment and Site effect) or variance of R<sub>sj</sub> if significant</td><td>0.11</td><td>0.51</td><td>0.12</td><td>0.24</td><td>0.41</td></tr><tr><td>Variance of ε<sub>kjj</sub> (Residuals)</td><td>0.004</td><td>0.20</td><td>0.24</td><td>0.35</td><td>0.50</td></tr><tr><td>R<sup>2</sup></td><td>0.53 (±0.04)</td><td>0.55 (±0.04)</td><td>0.55 (±0.04)</td><td>0.51 (±0.04)</td><td>0.47 (±0.05)</td></tr><tr><td>RMSE</td><td>0.06 (±0.003)</td><td>0.44 (±0.02)</td><td>0.50 (±0.02)</td><td>0.61 (±0.04)</td><td>0.75 (±0.05)</td></tr></tbody></table>

Model terms are given as bootstrapped mean and standard deviation (in brackets). Dataset used is indicated with the index: A when aborted experiments were removed, B, C, D when the $K_{sat}$ of aborted experiments were fixed at 10$^{-4}$, 10$^{-5}$, 10$^{-6}$ mm s$^{-1}$ respectively

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Distribution of the measured soil water content at saturation at 0–10 cm as a function of treatment and site</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Measured versus predicted values of the T-treatment with SITE\_mod model for $\theta_{sat}$ and $\log_{10}(K_{sat})$ in **A** aborted experiment exclusion case and **B** aborted experiment inclusion case with $K_{sat}$ fixed at 0.0001 mm s$^{-1}$ 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</figcaption></figure>

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.

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_{10}(K_{sat})$ between 0.50 and 0.75 when $K_{sat}$ of aborted tests was changed from $K_{sat}=10^{-4}$ to $K_{sat}=10^{-6}$ mm s$^{-1}$ (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^2$, i.e. $K_{sat}=10^{-4}$ mm s$^{-1}$ (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_{10}(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_{10}(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 $\theta_{sat}$ and $\log_{10}(K_{sat})$: rut depth, rut depth class ($\le 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 $\theta_{sat}$ and $\log_{10}(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 $\theta_{sat}$ and $\log_{10}(K_{sat})$ since the other variables proved to be less sensitive to soil compaction.

For PROX\_mod, both depth and intensity of skid degradation had an impact on $\theta_{sat}$ and $\log_{10}(K_{sat})$. In addition, their interactions depth\*intensity of skid trail degradation were significant for both $\theta_{sat}$ and $\log_{10}(K_{sat})$. The effect on $\theta_{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, $\theta_{sat}$ decreases when the soil degradation score increases. For $\log_{10}(K_{sat})$, the effect decreases when the soil degradation score increases at 0–10 cm. At 15–25 cm,
<table><caption>**Table 5** Relationship between the intensity of skid trail degradation, depth and soil hydraulic parameters while considering the site as random effect (= PROX_mod)</caption><thead><tr><th></th><th>θ<sub>sat</sub><sup>A</sup><br>m<sup>3</sup> m<sup>-3</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>A</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th><th>log<sub>10</sub>(K<sub>sat</sub>)<sup>B</sup><br>K<sub>sat</sub> in mm s<sup>-1</sup></th></tr></thead><tbody><tr><td>μ (Control at 0–10 cm)</td><td>0.64 (0.007)</td><td>-1.73 (0.05)</td><td>-1.72 (0.05)</td></tr><tr><td>α<sub>i</sub> (15–25 cm)</td><td>-0.06 (0.009)</td><td>-0.43 (0.07)</td><td>-0.44 (0.07)</td></tr><tr><td>γ<sub>i=1</sub> (Score degradation coefficient of level 1)</td><td>-0.023 (0.04)</td><td>-0.81 (0.21)</td><td>-1.21 (0.20)</td></tr><tr><td>γ<sub>i=2</sub> (Score degradation coefficient of level 2)</td><td>-0.044 (0.01)</td><td>-0.34 (0.10)</td><td>-0.48 (0.13)</td></tr><tr><td>γ<sub>i=3</sub> (Score degradation coefficient of level 3)</td><td>-0.010 (0.02)</td><td>-0.26 (0.11)</td><td>-0.26 (0.12)</td></tr><tr><td>α<sub>i</sub>;γ<sub>i=1</sub></td><td>0.005 (0.05)</td><td>0.64 (0.19)</td><td>0.50 (0.26)</td></tr><tr><td>α<sub>i</sub>;γ<sub>i=2</sub></td><td>-0.004 (0.01)</td><td>-0.16 (0.13)</td><td>-0.19 (0.16)</td></tr><tr><td>α<sub>i</sub>;γ<sub>i=3</sub></td><td>-0.070 (0.02)</td><td>-0.16 (0.14)</td><td>-0.16 (0.13)</td></tr><tr><td>Variance of R<sub>s</sub> (Site effect)</td><td>0.002</td><td>0.13</td><td>0.16</td></tr><tr><td>Variance of ε<sub>kij</sub> (Residuals)</td><td>0.004</td><td>0.20</td><td>0.25</td></tr><tr><td>R<sup>2</sup></td><td>0.53 (0.04)</td><td>0.59 (0.04)</td><td>0.59 (0.04)</td></tr><tr><td>RMSE</td><td>0.06 (0.003)</td><td>0.42 (0.02)</td><td>0.48 (0.03)</td></tr></tbody></table>

Model terms are given as bootstrapped mean and standard deviation (in brackets). Dataset used is indicated with the index: A when aborted experiments were removed (aborted experiment exclusion case), B when the K<sub>sat</sub> of aborted experiments were fixed at 10<sup>-4</sup> mm s<sup>-1</sup> (aborted experiment inclusion case)

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Measured versus predicted values of the T-treatment with PROX_mod model for θ<sub>sat</sub> and log<sub>10</sub>(K<sub>sat</sub>) in A aborted experiment exclusion case and B aborted experiment inclusion case with K<sub>sat</sub> fixed at 0.0001 mm s<sup>-1</sup> 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</figcaption></figure>

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<sub>sat</sub> (Tables 4 and 5, Figs. 7 and 8).

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

Predicted θ<sub>sat</sub> 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<sub>sat</sub> are mainly less than the reference values of K<sub>sat</sub> 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.

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 REFBEST scenario and to the observed data (Fig. 10). For
<figure><img src="image_13.png" /><figcaption>**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 log10(Ksat) of the T-treatment in A aborted experiment exclusion case and B aborted experiment inclusion case with Ksat fixed at 0.0001 mm s⁻¹ for aborted experiment</figcaption></figure>

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.

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

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 log10(Ksat). We expected that the scale parameter hg 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
<figure><img src="image_14.png" /><figcaption>**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</figcaption></figure>

<table><caption>**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</caption><thead><tr><th></th><th>REF_BEST</th><th>SITE</th><th>PROX</th><th>PTF_SITE</th><th>PTF_PROX</th></tr></thead><tbody><tr><td>10–30 cm</td><td>0.12</td><td>0.05</td><td>0.05</td><td>0.04</td><td>0.07</td></tr><tr><td>30–100 cm</td><td>0.20</td><td>0.10</td><td>0.09</td><td>0.09</td><td>0.07</td></tr></tbody></table>

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 ($\theta_{initial}$, $\theta_{sat}$) and all parameters estimated previously by the BEST method ($n$, $\eta$, $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*<sub>g</sub>. 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*<sub>sat</sub> for aborted tests slightly improved the assessment of the effect of compaction on *K*<sub>sat</sub>. 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*<sub>sat</sub> 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 $\theta_{sat}$ are calculated in the same way using soil bulk density and solid (or mineral) density. The low impact of traffic on $\theta_{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 redoxi-morphic features) and a high impact on $\theta_{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 $\theta_{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.

The use a visual score of skid trail degradation intensity did not improve the estimation of $\theta_{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 $\theta_{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 $\theta_{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.
<figure><img src="image_15.png" /><figcaption>**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</figcaption></figure>

**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*<sub>g</sub>, *n* and *η* were derived from the undisturbed soils (estimated with the BEST method or a PTF) while *K*<sub>sat</sub> and *θ*<sub>sat</sub> 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, Ksat 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 Ksat 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.

### Acknowledgements

This project would not have been possible without the technical support of our colleagues X. Montagny, C. Perinot, J. Levillain, J. Fiquepron and O. Monnier and interns A. Tornambe, M. Porentru and A. de Heaulme.

### Authors' contributions

Conceptualization: Manon Martin, Noémie Pousse, Stéphane Ruy; Methodology: Laurent Lassabatère; Formal analysis and investigation: Manon Martin, Noémie Pousse, Stéphane Ruy, André Chanzy, Laurent Lassabatère, Arnaud Legout; Writing – original draft preparation: Manon Martin; Writing – review and editing: Noémie Pousse, Stéphane Ruy, André Chanzy, Laurent Lassabatère, Arnaud Legout; Funding acquisition: Stephane Ruy, André Chanzy, Noémie Pousse; Supervision: Stephane Ruy, André Chanzy, Noémie Pousse.

### Funding

The experimental measurements were carried out within the EFFORTE project (Efficient forestry for sustainable and cost-competitive bio-based industry) which was financed by the Bio Based Industries Joint Undertaking under the European Union's Horizon 2020 research and innovation framework programme, grant agreement No 720712 and COPACEL. The experimental measurement of three additional sites ("ARF1", "ARF2" and "CAM") were funded by the French National Office of Forestry (Office National des Forêts, ONF). The two long-term monitored experimental sites "AZ25" and "CA" belong to French national research infrastructure ANAEE-France; we gratefully acknowledge the financial support successively received since 2007 from the French Ministry of Ecology and Sustainable Development (MEDD - Gessol), the French Ministry of Agriculture and Fisheries (MAP), GIP ECOFOR (AllEnvi-Soere), ANAEE-France (ANR-11-INBS-0001AnaEE-Services) and ONF. The BEERKAN measurements of the two long-term monitored experimental sites "AZ25" and "CA" and the statistical analysis of the dataset were carried out within the « GRAINE 2019 – VSoilForOAD», project which was funded by the French Agency for Ecological Transition (ADEME) under grant agreement n°2003 C 0060.

### Data availability

Data are described in a companion data-paper by Martin et al. (2024). Soil data used in this study are freely available at https://doi.org/10.15454/HS8U8D (Martin 2022).

## Declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

All the authors included in the manuscript agree to the publication of the article.

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹UMR EMMAH, INRAE, Avignon Université, 84000 Avignon, France. ²UMR LEHNA, Université de Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, 69518 Vaulx-en-Velin, France. ³UR BEF, INRAE, 54280 Champenoux, France.

⁴Département Recherche Et Développement, ONF, 73000 Chambéry, France.

Received: 10 January 2024 Accepted: 3 October 2024

Published: 7 November 2024

## References

Al Majou H, Bruand A, Duval O et al (2008) Prediction of soil water retention properties after stratification by combining texture, bulk density and the type of horizon. Soil Use Manag 24:383-391. https://doi.org/10.1111/j.1475-2743.2008.00180.x

Ampoorter E, Goris R, Cornelis WM, Verheyen K (2007) Impact of mechanized logging on compaction status of sandy forest soils. For Ecol Manage 241:162-174. https://doi.org/10.1016/j.foreco.2007.01.019

Angulo-Jaramillo R, Bagarello V, Iovino M, Lassabatère L (2016) Infiltration measurements for soil hydraulic Characterization. Springer

Angulo-Jaramillo R, Bagarello V, Di Prima S et al (2019) Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales. J Hydrol 576:239-261. https://doi.org/10.1016/j.jhydrol.2019.06.007

Assouline S (2006) Modeling the Relationship between Soil Bulk Density and the Hydraulic Conductivity Function. Vadose Zone Journal 5:697-705. https://doi.org/10.2136/vzj2005.0084

Aust WM, Burger JA, Carter EA et al (1998) Visually determined soil disturbance classes used as indices of forest harvesting disturbance. South J Appl for 22:245-250. https://doi.org/10.1093/sjaf/22.4.245

Bagarello V, Castellini M, Di Prima S, Iovino M (2014) Soil hydraulic properties determined by infiltration experiments and different heights of water pouring. Geoderma 213:492-501. https://doi.org/10.1016/j.geoderma.2013.08.032

Ballard TM (2000) Impacts of forest management on northern forest soils. For Ecol Manage 133:37-42. https://doi.org/10.1016/S0378-1127(99)00296-0

Bates D, Mächler M, Bolker B, Walker S (2015) Fitting Linear Mixed-Effects Models Using lme4. J Stat Softw 67:1-48. https://doi.org/10.18637/jss.v067.i01

Bottinelli N, Hallaire V, Goutal N et al (2014) Impact of heavy traffic on soil macroporosity of two silty forest soils: Initial effect and short-term recovery. Geoderma 217-218:10-17. https://doi.org/10.1016/j.geoderma.2013.10.025

Braud I, De Condappa D, Soria JM et al (2005) Use of scaled forms of the infiltration equation for the estimation of unsaturated soil hydraulic properties (the Beerkan method). Eur J Soil Sci 56:361-374. https://doi.org/10.1111/j.1365-2389.2004.00660.x

Brooks R, Corey C (1964) Hydraulic properties of porous media. Fort Collins, Colorado state

Burdine NT (1953) Relative Permeability Calculations From Pore Size Distribution Data. J Petrol Technol 5:71-78. https://doi.org/10.2118/225-G

Dexter AR (2004) Soil physical quality: Part I. Theory, effects of soil texture, density, and organic matter, and effects on root growth. Geoderma 120:201-214. https://doi.org/10.1016/j.geoderma.2003.09.004

Di Prima S, Concialdi P, Lassabatère L et al (2018) Laboratory testing of Beerkan infiltration experiments for assessing the role of soil sealing on water infiltration. CATENA 167:373-384. https://doi.org/10.1016/j.catena.2018.05.013
Di Prima S, Stewart RD, Abou Najm MR et al (2021) BEST-WR: An adapted algorithm for the hydraulic characterization of hydrophilic and water-repellent soils. J Hydrol 603:126936. https://doi.org/10.1016/j.jhydrol.2021.126936

Driscoll CT, Schecher WD (1990) The chemistry of aluminum in the environment. Environ Geochem Health 12:28-49. https://doi.org/10.1007/BF01734046

Ebeling C, Fründ H-C, Lang F, Gaertig T (2017) Evidence for increased P availability on wheel tracks 10 to 40years after forest machinery traffic. Geoderma 297:61-69. https://doi.org/10.1016/j.geoderma.2017.03.003

Efron B, Tibshirani RJ (1993) An introduction to the bootstrap. In: Monographs on statistics and applied probability, Springer Sciences-business media edition

Findeling A, Garnier P, Coppens F et al (2007) Modelling water, carbon and nitrogen dynamics in soil covered with decomposing mulch. Eur J Soil Sci 58:196-206. https://doi.org/10.1111/j.1365-2389.2006.00826.x

Frey B, Kremer J, Rüdt A et al (2009) Compaction of forest soils with heavy logging machinery affects soil bacterial community structure. Eur J Soil Biol 45:312-320. https://doi.org/10.1016/j.ejsobi.2009.05.006

Fründ H-C, Averdick A (2016) Soil aeration and soil water tension in skidding trails during three years after trafficking. For Ecol Manage 380:224-231. https://doi.org/10.1016/j.foreco.2016.09.008

Goutal N, Boivin P, Ranger J (2012) Assessment of the natural recovery rate of soil specific volume following forest soil compaction. Soil Sci Soc Am J 76:1426-1435. https://doi.org/10.2136/sssaj2011.0402

Goutal N, Keller T, Défossez P, Ranger J (2013) Soil compaction due to heavy forest traffic: Measurements and simulations using an analytical soil compaction model. Ann for Sci 70:545-556. https://doi.org/10.1007/s13595-013-0276-x

Goutal-Pousse N, Lamy F, Ranger J, Boivin P (2016) Structural damage and recovery determined by the colloidal constituents in two forest soils compacted by heavy traffic. Eur J Soil Sci 67:160-172. https://doi.org/10.1111/ejss.12323

Greacen E, Sands R (1980) Compaction of forest soils. A Review. Australian J Soil Res 18:163-189 https://doi.org/10.1071/SR9800163

Hartshorne R, Dicken SN (1935) A classification of the agricultural regions of Europe and North America on a uniform statistical basis. Ann Assoc Am Geogr 25:99-120. https://doi.org/10.2307/2560652

Haverkamp R, Ross PJ, Smettem KRJ, Parlange JY (1994) Three-dimensional analysis of infiltration from the disc infiltrometer: 2. Physically-Based Infiltration Equation Water Resources Research 30:2931-2935. https://doi.org/10.1029/94WR01788

Haverkamp R, Parlange J-Y, Starr JL, et al (1990) Infiltration under ponded conditions. 3. A predictive equation based on physical parameters. Soil Sci 149:292. https://journals.lww.com/soilsci/Abstract/1990/05000/INFILTRATION_Under_ponded CONDITIONS__3__A.6.aspx

Heninger R, Scott W, Dobkowski A et al (2002) Soil disturbance and 10-year growth response of coast Douglas-fir on nontilled and tilled skid trails in the Oregon Cascades. Can J for Res 32:233-246. https://doi.org/10.1139/x01-195

Keller T, Sandin M, Colombi T et al (2019) Historical increase in agricultural machinery weights enhanced soil stress levels and adversely affected soil functioning. Soil Tillage Res 194:1-12. https://doi.org/10.1016/j.still.2019.104293

Kozlowski TT (1999) Soil Compaction and Growth of Woody Plants. Scand J for Res 14:596-619. https://doi.org/10.1080/02827589908540825

Lafolie F (1991) Modelling water flow, nitrogen transport and root uptake including physical non-equilibrium and optimization of the root water potential. Fertilizer Res 27:215-231. https://doi.org/10.1007/BF01051129

Lassabatere L, Angulo-Jaramillo R, Soria Ugalde JM et al (2006) Beerkan Estimation of Soil Transfer Parameters through Infiltration Experiments-BEST. Soil Sci Soc Am J 70:521-532. https://doi.org/10.2136/sssaj2005.0026

Lassabatere L, Angulo-Jaramillo R, Soria-Ugalde JM et al (2009) Numerical evaluation of a set of analytical infiltration equations. Water Resour Res 45:W12415. https://doi.org/10.1029/2009WR007941

Lassabatere L, Di Prima S, Angulo-Jaramillo R et al (2019) Beerkan multi-runs for characterizing water infiltration and spatial variability of soil hydraulic properties across scales. Hydrol Sci J 64:165-178. https://doi.org/10.1080/02626667.2018.1560448

Lassabatere L, Peyneau P-E, Yilmaz D et al (2021) A scaling procedure for straightforward computation of sorptivity. Hydrol Earth Syst Sci 25:5083-5104. https://doi.org/10.5194/hess-25-5083-2021

Lassabatere L, Peyneau P-E, Yilmaz D et al (2023) Mixed formulation for an easy and robust numerical computation of sorptivity. Hydrol Earth Syst Sci 27:895-915. https://doi.org/10.5194/hess-27-895-2023

Mariotti B, Hoshika Y, Cambi M et al (2020) Vehicle-induced compaction of forest soil affects plant morphological and physiological attributes: A meta-analysis. For Ecol Manage 462:118004. https://doi.org/10.1016/j.foreco.2020.118004

Martin M, Chanzy A, Lassabatere L, et al (2024) Hydraulic properties for a wide range of undisturbed and compacted French forest soils: in situ measurements and estimation with the BEST method. Ann Forest Sci. https://doi.org/10.1186/s13595-024-01262-7

Martin M (2019) Développement d'un modèle de transfert hydrique des sols forestiers partiellement tassés dans un contexte de données parcimonieuses. Thèse de doctorat, Université d'Avignon, https://tel.archives-ouvertes.fr/tel-02518560

Martin M (2022) A hydraulic properties dataset to provide the water flow of compacted forest soils. V7. Recherche Data Gouv. https://doi.org/10.15454/HS8U8D

McNabb DH, Startsev AD, Nguyen H (2001) Soil Wetness and Traffic Level Effects on Bulk Density and Air-Filled Porosity of Compacted Boreal Forest Soils. Soil Sci Soc Am J 65:1238-1247. https://doi.org/10.2136/sssaj2001.6541238x

Mohieddinne H, Brasseur B, Spicher F et al (2019) Physical recovery of forest soil after compaction by heavy machines, revealed by penetration resistance over multiple decades. For Ecol Manage 449:117472. https://doi.org/10.1016/j.foreco.2019.117472

Naderi-Boldaji M, Keller T (2016) Degree of soil compactness is highly correlated with the soil physical quality index S. Soil Tillage Res 159:41-46. https://doi.org/10.1016/j.still.2016.01.010

Pachepsky YA, Rawls WJ (2003) Soil structure and pedotransfer functions. Eur J Soil Sci 54:443-451. https://doi.org/10.1046/j.1365-2389.2003.00485.x

Pousse N, Bonnaud P, Legout A et al (2022) Forest soil penetration resistance following heavy traffic: A 10-year field study. Soil Use Manag 38:815-835. https://doi.org/10.1111/sum.12730

Ranger J, Bonnaud P, Santenoise P et al (2021) Effect of limited compaction on soil solution chemistry in two acidic forest ecosystems: Changes, recovery and impact of liming. For Ecol Manage 499:119538. https://doi.org/10.1016/j.foreco.2021.119538

Richard G, Cousin I, Sillon JF et al (2001) Effect of compaction on the porosity of a silty soil: influence on unsaturated hydraulic properties. Eur J Soil Sci 52:49-58. https://doi.org/10.1046/j.1365-2389.2001.00357.x

Richards LA (1931) CAPILLARY CONDUCTION OF LIQUIDS THROUGH POROUS MEDIUMS. Physics 1:318-333. https://doi.org/10.1063/1.1745010

Salmivaara A, Launiainen S, Perttunen J et al (2021) Towards dynamic forest trafficability prediction using open spatial data, hydrological modelling and sensor technology. Forestry 93:662-674. https://doi.org/10.1093/FORESTRY/CPAA010

Schaap MG, Leij FJ, van Genuchten MTh (2001) rosetta : a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J Hydrol 251:163-176. https://doi.org/10.1016/S0022-1694(01)00466-8

Schack-Kirchner H, Fenner PT, Hildebrand EE (2007) Different responses in bulk density and saturated hydraulic conductivity to soil deformation by logging machinery on a Ferralsol under native forest. Soil Use Manag 23:286-293. https://doi.org/10.1111/j.1475-2743.2007.00096.x

Smettem KRJ, Parlange JY, Ross PJ, Haverkamp R (1994) Three-dimensional analysis of infiltration from the disc infiltrometer: 1. A Capillary-Based Theory Water Resources Research 30:2925-2929. https://doi.org/10.1029/94WR01787

Smith R, Ellies A, Horn R (2000) Modified Boussinesq's equations for non uniform tire loading. J Terrramech 37:207-222. https://doi.org/10.1016/S0022-4898(00)00007-0

Tóth B, Weynants M, Nemes A et al (2015) New generation of hydraulic pedotransfer functions for Europe: New hydraulic pedotransfer functions for Europe. Eur J Soil Sci 66:226-238. https://doi.org/10.1111/ejss.12192

van Genuchten MTh (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 44:892-898. https://doi.org/10.2136/sssaj1980.03615995004400050002x
Vega-Nieva DJ, Murphy PNC, Castonguay M et al (2009) A modular terrain model for daily variations in machine-specific forest soil trafficability. Can J Soil Sci 89:93-109. https://doi.org/10.4141/CJSS06033

Wessolek G, Kaupenjohann M, Renger M (2009) Bodenphysikalische Kennwerte und Berechnungsverfahren für die Praxis. Technische Universität Berlin, Berlin. https://doi.org/10.13140/RG.2.2.27053.10729

Weynants M, Montanarella L, Toth G, et al (2013) European HYdropedological Data Inventory (EU-HYDI), GD Luxembourg: Publications Office of the European Union. Luxembourg, https://data.europa.eu/doi/https://doi.org/10.2788/5936

Williamson JR, Neilsen WA (2000) The influence of forest site on rate and extent of soil compaction and profile disturbance of skid trails during ground-based harvesting. Can J for Res 30:1196-1205. https://doi.org/10.1139/x00-041

Wösten JHM, Lilly A, Nemes A, Le Bas C (1999) Development and use of a database of hydraulic properties of European soils. Geoderma 90:169-185. https://doi.org/10.1016/S0016-7061(98)00132-3

Wösten JHM, Pachepsky YA, Rawls WJ (2001) Pedotransfer function: bridging the gap between available basic soil data and missing soil hydraulic characterics. J Hydrol 123-150. https://doi.org/10.1016/S0022-1694(01)00464-4

Yilmaz D, Lassabatere L, Angulo-Jaramillo R et al (2010) Hydrodynamic Characterization of Basic Oxygen Furnace Slag through an Adapted BEST Method. Vadose Zone J 9:107. https://doi.org/10.2136/vzj2009.0039

Yilmaz D, Di Prima S, Stewart RD et al (2022) Three-term formulation to describe infiltration in water-repellent soils. Geoderma 427:116127. https://doi.org/10.1016/j.geoderma.2022.116127

Ziesak M (2003) Avoiding soil damages, caused by forest machines. In Wide, M.I. and Hallberg, I., Eds. Proceedings: 2nd Forest Engineering Conference, 12-15 May 2003, Vaxjö, Sweden, Skogforsk Arbetsrapport No 536:3-11

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.