Hydraulic properties for a wide range of undisturbed and compacted French forest soils: in situ measurements and estimation with the BEST method

Abstract

Key message

The dataset provides hydraulic properties estimated using the Beerkan Estimation of Soil Transfer (BEST) method, on undisturbed and on compacted and rutted French forest soils. It allows a reliable assessment of the effect of traffic on soil permeability. However, hydraulic properties could not be estimated on extremely rutted soils, underscoring the necessity for tailored protocols for these conditions.

1 Background

Traffic of forestry machines generates, from the first passage on, soil degradation by compaction and rutting (McNabb et al. 2001). This is a long-term degradation as natural recovery is very low in forest soils (Mohieddinne et al. 2019). To avoid circulation over the entire plot, skid trails (i.e., permanent corridors dedicated to the circulation of forestry machinery) are now routinely set up in French forests. However, skid trails must remain practicable in the long term, and forest managers must choose the logging sites and stop traffic according to the weather conditions.

The soil bearing capacity depends on its composition (texture, bulk density, organic matter content, structure) and its moisture (McNabb et al. 2001), which depends on the hydraulic properties and climatic conditions. Mechanistic soil water flow models can predict soil water dynamics and may be integrated in decision support tools for planning forestry interventions. Although soil hydraulic properties are key parameters of these models and are good indicators of soil health (https://environment.ec.europa.eu/topics/soil-and-land/soil-health_en), few data exist in a forest context and even less on compacted soils in forest ecosystems. Indeed, published datasets such as USDA (Hartshorne and Dicken 1935) and HYPRES (Wösten et al. 1999) cover a small fraction of forest soils without any data for rutted soils.

The dataset built up in this study is composed of 417 Beerkan infiltration experiments over 19 different plots with measurements of bulk density, initial and final water contents, and particle size distributions for undisturbed soils (control, C-treatment) and soils underneath skid trails (trafficked, T-treatment). The hydraulic properties were determined with the Beerkan Estimation of Soil Transfer parameter (BEST) methods that are reviewed in Angulo-Jaramillo et al. (2019). BEST uses the van Genuchten formulation (van Genuchten 1980) with the Burdine condition (Burdine 1953) to describe the water retention curve and the Brooks and Corey formulation (Brooks and Corey 1964) to describe the hydraulic conductivity curve. The scale parameters (h_g, K_sat in Eqs. 1 and 2) are estimated from fitting the analytical models developed by Haverkamp et al. (1994) to the observed cumulative infiltration with five different strategies corresponding to the five BEST methods. The shape parameters (n, m, and η in Eqs. 1 and 2) are estimated from the particle size distribution and pedotransfer functions specific to the BEST methods (Lassabatere et al. 2006).

2 Methods

2.1 Description of experimental plots

Nineteen plots in 16 forests mainly located in North-Eastern France were studied. Data were collected on plain forest stands where skid trails had already been installed. Site location and related general characteristics are presented in Table 1. The skid trails have been circulated at least once, and the next harvest was planned for the next 2 years after the water infiltration measurements. The forests have been chosen to cover a large range of soil textures with a majority of silty loam texture, which is known to be highly sensitive to soil compaction (Fig. 1).

Table 1 Description of the experimental plots studied: name, location, stand composition, and soil textural class according to the USDA classification (Hartshorne and Dicken 1935). FD Forêt Domaniale (state_owned forest), SiCl silty clay, Cl clay, SaClLo sandy clay loam, SaLo sandy loam, SiLo silty loam, SiClLo silty clay loam

Fig. 1

Texture of the undisturbed soils studied according to the USDA classification (US Department of Agriculture classification system) and European classification (Hartshorne and Dicken 1935). The arrows along the axes of the soil texture triangle indicate the directions to assess the percentages of sand, silt, and clay

2.2 Soil sampling and infiltration experiments

Between 4 and 8 infiltration measurements were performed, after removing the forest floor, in the 0–10 cm and the 15–25 cm soil layers for each treatment and each plot, depending on the number of unsuccessful infiltration tests. Those corresponded to tests where the time required for infiltration was excessively long, preventing measurement of at least five cumulative infiltration times within 4 h. We carried out as many infiltration tests as needed to get at least three successful experiments per soil layer, site, and treatment. Measurements were also done in the 30–40 cm soil layer on three plots, but a lot of them failed due to very long infiltration time, especially in the T-treatment. Therefore, these measurements are considered as unsuccessful and thus not presented. Measurements were performed in the center of the wheel tracks of the skid trails (T-treatment) and in the undisturbed forest area close to the skid trails (C-treatment). Based on weather conditions prior to the experimental campaigns, infiltration measurements were performed when the soil was not too dry to avoid water repellency and when it was not too wet to measure both transient and steady states (Lassabatere et al. 2009). Indeed, when the soil is initially close to saturation, steady state is quickly reached and the transient state is poorly described. Prior to the Beerkan infiltration experiments, two soil core samples were extracted to record the bulk density and the initial water content. A measure of the final water content was also carried out at the end of the Beerkan infiltration experiments. Those consist of recording the infiltration time of successive fixed volumes of water (here between 2.5 and 8 mm of water depths) poured into a simple annular ring between 10 and 19.6 cm in diameter. According to the measurement campaigns, different cylinders were used. The smallest cylinders and the highest water volumes used for the first measurement campaigns lead to cumulative infiltrations with less than a dozen points, which is not enough to ensure the proper use of the BEST methods. Therefore, larger cylinders and smaller water volumes were used to increase the BEST method’s reliability for the following experiments. The experiment is considered as completed when 3 successive volumes of water lead to similar infiltration times, meaning that steady state is reached. In addition to infiltration experiments, a composite sample has been made per treatment and layer, with soil samples taken next to the infiltration experiments to determine particle size distribution, pH, organic matter, and carbonate content.

2.3 Hydraulic characterization

The hydraulic properties were determined with the BEST method (Lassabatere et al. 2006). BEST uses the van Genuchten formulation (van Genuchten 1980) with the Burdine condition (Burdine 1953) to describe the water retention curve (Eq. 1) and the Brooks and Corey formulation (Brooks and Corey 1964) to describe the hydraulic conductivity curve (Eq. 2):

$$\frac{\theta -{\theta }_{r}}{{\theta }_{sat}-{\theta }_{r}}={\left(1+{\left(\frac{h}{{h}_{g}}\right)}^{n}\right)}^{-m} \text{with}\ m=1-{}^{2}\!\left/ \!{}_{n}\right.$$

(1)

$$\frac{K(\theta )}{{K}_{sat}}={\left(\frac{\theta -{\theta }_{r}}{{\theta }_{sat}-{\theta }_{r}}\right)}^{\eta } \text{with }\eta =\frac{2}{nm}+2+p$$

(2)

where θ_sat and θ_r are respectively the soil volumetric water content at saturation and the residual volumetric water content, h_g is a scale parameter for water pressure head (mm), K_sat is the hydraulic conductivity at saturation (mm s⁻¹), p is a tortuosity parameter taken as p = 1 for Burdine condition, and n, m, and η are shape parameters.

BEST algorithms estimate the shape parameters from the particle size distribution considering the pedotransfer functions described in Lassabatere et al. (2006). θ_sat is deduced from the bulk density (ρ_b) measurement and the solid density (ρ_s = 2.65 g cm⁻³): θ_sat = 1 − ρ_b / ρ_s. Finally, θ_r is assumed to be null (Lassabatere et al. 2006). The scale parameters K_sat and h_g are estimated from the inversion of the cumulative infiltration curves I(t) (mm) by fitting the analytical models developed by Haverkamp et al. (1994) to the observations using five different BEST algorithms. BEST-slope, BEST-intercept, and BEST-steady algorithms are based on the same mathematical framework but their fitting process differs; BEST-WR and BEST-WR-3T are based on a modified mathematical framework of BEST-slope to account for water repellency (Abou Najm et al. 2021). For the sake of simplicity and to avoid numerical problems, all BEST algorithms rely on the approximate expansions developed by Haverkamp et al. (1994) for transient and steady states instead of the quasi-exact implicit model developed by the same authors (Lassabatere et al. 2009). The different BEST algorithms might be summarized as follows:

• BEST-slope fits the analytical model to the experimental data using the slope of the linear regression of the steady state part of I(t) as a constraint between sorptivity and hydraulic conductivity (Lassabatère et al. 2006): K_sat = q_+∞ − A S²

• BEST-intercept fits the model to the experimental data using the intercept of the linear regression of the steady state part of I(t) as a constraint between sorptivity and hydraulic conductivity (Yilmaz et al. 2010): K_sat = C S²/b_+∞

• BEST-steady fits the model to the experimental data using both the slope and the intercept of the linear regression of the steady state part of I(t) and making no uses of the transient part (Bagarello et al. 2014): K_sat = q_+∞ − A S² and K_sat = C S²/b_+∞

• BEST-WR relies on the correction of water infiltration rate proposed by Abou Najm et al. (2021) for water repellent soils. Abou Najm et al. (2021) combined their correction factor with the two-term approximate expansion proposed by Haverkamp et al. (1994) for transient state. Di Prima et al. (2021) were the first authors to implement the Abou Najm’s corrected expansion with the BEST approach to build up the BEST-WR method dedicated to water repellent soils (Di Prima et al. 2021). Water repellence is characterized by a new scaling factor, α_WR. Afterwards, the other hydraulic parameters are estimated as in the BEST-slope method. The BEST-WR algorithm leads to better fit for both hydrophilic (concave cumulative infiltration curve) and hydrophobic (convex cumulative infiltration curve) soils.

• BEST-WR-3T extends the BEST-WR method to increase its accuracy for Beerkan tests for long tests performed on highly water repellent soils (Yilmaz et al. 2022). In more details, Yilmaz et al. (2022) combined the correction factor proposed by Abou Najm et al. (2021) with the more precise three-term approximate expansions developed by Haverkamp et al. (1994) for transient state, whereas the BEST-WR considered only the two first terms. This new method proved more accurate under given circumstances.

Angulo-Jaramillo et al. (2019) recommended to combine the three algorithms BEST-slope, BEST-intercept, and BEST-steady to overcome fitting errors and to estimate accurate h_g and K_sat parameters. We extended this recommendation to BEST-WR and BEST-WR-3T algorithms as their fitting errors proved lower than BEST-slope, BEST-intercept, and BEST-steady (Fig. 2), especially for hydrophobic soils (Fig. 3, examples of improvements with the BEST-WR methods). Besides, BEST-WR and BEST-WR-3T estimate similar values for S and K_sat in most cases, and, for severely water-repellent soils, only BEST-WR-3T yields robust fits to experimental data. Therefore, we performed the five algorithms by using encoded BEST (Lassabatere et al. 2013) on Scilab open source software (Campbell et al. 2006). The quality of the fit for the 5 algorithms was assessed using the relative error Er (Eq. 3) and the bias (Eq. 4). To ensure accurate modeling of experimental data, for each infiltration experiment, we kept only the fits that complied with two thresholds at the same time, i.e., for the relative error and the bias (Lassabatère et al. 2006):

$$Er=100*\sqrt{\frac{{\sum }_{i=1}^{k}{\left({I}_{exp,i}-{I}_{sim,i}\right)}^{2}}{{\sum }_{i=1}^{k}{{I}_{exp,i}}^{2}}}$$

(3)

$$Bias=\frac{1}{k}{\sum }_{i=1}^{k}\left({I}_{sim,i}-{I}_{exp,i}\right)$$

(4)

where I_exp,i and I_sim,i are, respectively, observed and simulated cumulative infiltrations at each measurement point i of a given infiltration run, and k is the number of measurement points per infiltration experiment. We also controlled that the fitted model accurately reproduced the shape of the experimental cumulative infiltrations, in particular its concavity or convexity.

Fig. 2

Relative error of fit (Eq. 3), for experimental data displaying relative errors below 10%, as a function of the BEST algorithm (intercept, slope, steady, WR, and WR-3T), the soil layer (0–10 and 15–25 cm), and the treatment (C-treatment being the undisturbed soil, T-treatment being the compacted soil under skid trails)

Fig. 3

Examples of cumulated water height infiltrated as a function of time: “ABB” plot at 0–10 cm depth, C-treatment (top) and T-treatment (bottom) observations as black dots, BEST simulations as blue line for two algorithms BEST-slope on the left and BEST-WR on the right

2.4 Dataset content

The dataset includes (i) an overview of the bare characteristics of the soil plot, i.e., particle size distribution, pH, organic matter, and carbonate content, (ii) cumulated infiltration curves measured on the experimental sites during the period 2017–2021, (iii) bulk densities, initial and final soil water contents measured respectively before and after the infiltration experiments, and (iv) hydraulic properties estimated with the five algorithms of the BEST method and the chosen criteria for goodness of fits. Each infiltration experiment is identified by a specific ID following this scheme: PLOT-TREATMENT-DEPTH-REPLICATE.

3 Access to the data and metadata description

The data and metadata description are available (Martin 2022): https://doiorg.publicaciones.saludcastillayleon.es/10.15454/HS8U8D.

4 Technical validation

4.1 Aborted infiltration tests

Some infiltration experiments were considered as failed since the time needed to infiltrate was too long and did not allow measuring at least five cumulative infiltration times in 4 h. They concerned 18 infiltration experiments over the 417 (4.3%). All are found in the skid trail (T-treatment), and most of them in the AZ site. Soil texture, organic matter content, or initial water content did not explain Beerkan measurement failures. They seemed to be due to very low water infiltration rates (maximum cumulated infiltration height divided by maximum infiltration time), ranging from 0.0003 to 0.003 mm s⁻¹, whereas the water infiltration rates were much higher ranging from 0.001 to 0.60 mm s⁻¹ for the remaining 399 infiltration tests. All were performed on skid trails and reflect a significant impact of compaction. These tests could not be treated with BEST methods and did not provide any estimate for soil hydraulic curves. The exclusion of these tests from the database used to assess the compaction effect on hydraulic properties of the skid trails may have introduced bias, potentially underestimating the impact of soil compaction (Martin et al. 2024). Indeed, the discarded tests corresponded to an extreme or at least very important effect of compaction. All these tests are referred to as “aborted tests” in the following.

4.2 Selection of infiltration experiments based on relative errors and biases

Our aim was to keep the least biased possible experimental dataset to further study the effect of traffic (Martin et al. 2024), i.e., to keep only the BEST fits displaying no systematic over- or underfitting.

For the remaining 399 infiltration experiments, fitting BEST algorithms to experimental data lead to relative errors ranging between 0.2 and 446%, and bias ranging between − 83 and 152 mm. All the fits with relative errors Er smaller than or equal to 5% were automatically kept (Lassabatère et al. 2006). These runs were associated to low bias values (Fig. 4). All the fits with relative error, Er, higher than or equal to 10% were automatically discarded because of the very poor quality of fits (Figs. 3 and 4) even if the associated bias may be low, similar to the ones obtained with a good quality of fit (Er less than 5%). Indeed, a given model may be good for the prediction of the average without properly modeling the shape of the curve. Between these two thresholds of relative errors, some fits displayed general shapes without systematic over- or underfitting (low bias). Besides, in the T-treatment, a higher percentage of BEST simulations displayed relative errors above 5% than in the C-treatment (Table 2). Removing all the BEST simulations with relative errors above 5% could lead to an underestimation of the effect of traffic on hydraulic properties since this selection simultaneously removes more tests from the trafficked area than the control area and removes the tests associated to the higher levels of compaction.

Fig. 4

Bias as a function of the relative error, for all the BEST fits displaying relative errors below 100% and bias comprised between − 50 and 50 mm, per BEST algorithm. The blue horizontal lines represent the minimum and maximum of bias for the infiltration tests displaying relative errors lower than or equal to 5%; the green and the red vertical line indicate, respectively, 5 and 10% relative error

Table 2 Percentages of successful BEST simulations per treatment (C-treatment being the undisturbed soil, T-treatment being the compacted soil under skid trails), soil layer (0–10 and 15–25 cm), and class of relative error (Eq. 3) and bias (Eq. 4) of the fit for the 5 algorithms

Consequently, we chose the minimum and maximum bias values of the fits displaying relative errors smaller or equal to 5%, i.e., − 2.91 mm and 2.57 mm, to select and maintain the errors of the fit between 5 and 10%. With these thresholds on bias and relative error, fitting experimental data with BEST algorithms lead to poor quality of fit with all the five methods for 84 infiltration tests (20%). Those were discarded and referred to as “rejected tests or simulations” in the following.

4.3 Field and simulation success rate

We could not find any significant relationship between success rates of Beerkan runs (field) and BEST simulations on the one hand, and soil layer, particle size distribution, organic matter content, bulk density, infiltration rate (total cumulated infiltrated water divided by cumulated time; Fig. 5), initial water content (Fig. 6), and treatment, on the other hand, even for the two plots where particle size distribution and organic matter content were measured at each infiltration measurement point (“AZ25” and “CA”). An assessment of soil structure might have improved the explanation of success rates. The most important variables explaining success rate were the infiltration rate (IR) and the initial water content. The higher the infiltration rate and the lower the initial water content, the higher the chance for success for both Beerkan runs and BEST simulations (Fig. 5). Yet the link was neither tight nor precise, as for low infiltration rates, success was as likely as failure for Beerkan runs (field) and BEST simulations (high simulations errors; Fig. 5). Aborted Beerkan runs and rejected BEST simulations are associated with high initial water content (Fig. 6) for 6 plots (ABB, BTG, FTG, FTG148, SAU, SRG), among which two have 100% of simulations failure (FTG both treatments and BTG T-treatment). However, this was not always the case: for example, the initial water content of the plot “H” with a high BEST simulation success rate is as high as the one of “BTG” with a low BEST simulation success rate. Therefore, a failure in the BEST simulation is not always due to a high initial water content, but may also be the result of a high level of concavity (see Fig. 7, FTG plot and BTG plot T-treatment) that cannot be modeled by any of the different BEST algorithms. So far, we addressed the case of convexity. When too strong, water repellence may either prevent from any Beerkan test (because of too small amount of infiltrated water and too long infiltration time) or prevent from treating the cumulative infiltration with BEST-WR algorithms because of the convexity of the observed cumulative infiltrations. The opposite situation may occur, when the shape of the experimental data shows an extreme concavity. Under these circumstances, the model underestimates the observed cumulative infiltration curve at short times. A large concavity might be caused by the formation of a thin sealed layer at the surface (Di Prima et al. 2018), air entrapment, or soil layering (Alagna et al. 2016). These conditions prevent from the use of the regular BEST methods. They are bound to soil structure and may be poorly predicted by particle size distribution, bulk density, and organic matter content. Even if the BEST-WR and the BEST-WR-3T models improve the goodness of fit for convex curves, these are not appropriate for very concave curves and concave-then-convex curves. All BEST methods including BEST-WR and BEST-WR-3T had some difficulties when the observations presented a high concavity at the beginning of the infiltration measurement (Yilmaz et al. 2022).

Fig. 5

Initial volumetric water content (n = 417) plotted against infiltration rate (IR, total cumulated infiltrated volume of water divided by the cumulated time), segmented by fit quality into distinct groups: measurement failure (< 5 water volumes in 4 h), successful fits, and poor fits indicating relatively high relative error and bias

Fig. 6

Initial volumetric water content (n = 417) per plot and success rate of water infiltration measurement with the Beerkan method (less than 5 water volumes in 4 h = failure of in situ measurements) and of simulation with the BEST methods (poor BEST fit = relative error and bias of the BEST fits too high)

Fig. 7

Infiltration curves (infiltrated water height as a function of time) measured immediately below the litter (0–10 cm) per site (site code in the label of each group of infiltration curves) and treatment with undisturbed soil (C-treatment) in green and soil under skid trails (T-treatment) in Burgundy

5 Reuse potential and limits

5.1 Potentialities and novelty of the dataset

Our dataset includes forest stand and soil characteristics as well as infiltration curves required to estimate the hydraulic parameters with the BEST method. Forest soils represent only a few percentages of international datasets as USDA (Hartshorne and Dicken 1935) or HYPRES (Wösten et al. 1999): our dataset could thus help to develop pedotransfer functions adapted to the specificities of forest soils.

In addition, it compares undisturbed with trafficked soils through a large range of textural classes, which is very original and was scarcely addressed in the literature (Martin et al. 2024). Moreover, the use of a single and simple method to determine the hydraulic properties allows the increase of the number of replicates, thus facilitating the reproducibility of the measurements and decreasing estimate uncertainty.

It can also be used from a more applied viewpoint by developing soil degradation indicators or decision-support tools regarding the planning of silvicultural interventions (trafficking, harvest, etc.) or drought prediction and water management.

Then, the water cumulative infiltrations of the proposed dataset can complement other databases like the Soil Water Infiltration Global database that collects infiltration curves from all over the world (Rahmati et al. 2018). Further investigations may address the comparison of our data with datasets of cumulative infiltrations obtained for forest sites.

5.2 Limits of the dataset

It is very uneasy to gather the relevant information about the history of logging for the different sites, which leads to uncertainty about the absence of traffic on control treatments.

The limits of the method used to estimate the hydraulic properties are reached when the permeability of the soil is very low, especially on compacted soils, leading to unusable infiltration curves. Similar problems rise when the soil is very water-repellent or self-sealing, conveying to the observed cumulative infiltrations convexity or extra concavity and spoiling the treatment with BEST methods. Consequently, the hydraulic parameters were obtained only for the tests performed on sites where the effect of compaction was not too extreme, which minimizes our conclusions on the impact of soil compaction on the soil hydraulic conductivity at saturation. Hydraulic parameters could not be estimated from these infiltration curves but they revealed very slow saturated hydraulic conductivity. Specific and appropriate devices should be deployed for the characterization of very low permeability soils.