Hydrodynamic Parameter Estimation for Simulating Soil-Vegetation-Atmosphere Hydrology Across Forest Stands in the Strengbach Catchment

1. Introduction

Bibliometric analyses reveal a growing interest in the study of forest ecosystems, likely reflecting the increasing challenges they face [1]. Indeed, on the one hand, forests provide a wide range of ecosystem services, including food production, timber and fuel supply, water conservation and regulation, carbon sequestration, biodiversity protection, climate regulation, as well as ecotourism and other recreation activities [2]. However, one the other hand, their long-term sustainability is increasingly threatened by multiple anthropogenic pressures, particularly the expansion of agricultural lands and urban areas, along with the adverse impacts of climate change (CC) [3]. CC, in turn, manifests through prolonged droughts, disruptions in water availability, soil erosion, and declining soil fertility - factors that collectively weaken forest ecosystems, reduce their capacity to withstand and recover from such perturbations [3,4] and ultimately increases global tree mortality [5]. A combination of interacting factors is at play; for instance, the more frequent droughts observed over the past decade have not only caused water shortages and induced physiological stress in trees, but have also triggered indirect damage through pest outbreaks - most notably bark beetle infestations [6]- and increased the risk of large-scale fires [7].

Water is a key component to understand ecosystem functioning. Its circulation and the complex interactions within the soil–tree–atmosphere continuum are central concerns in ecohydrology. While climate and soil moisture regulate vegetation dynamics, forest ecosystems also strongly influence the overall water balance and generate multiple feedback loops [8]. Hence, soil moisture and vapor pressure deficit appear among the most influential factors controlling forest water use [9].

Characterizing the spatial and temporal dynamics of soil water storage is therefore essential. This need has, for instance, motivated the development of large-scale soil moisture monitoring databases [10] and efficient treatment of satellite-based measurements [11]. Ground-based observation networks, focusing on the critical zone and supporting interdisciplinary research, continue to play a crucial role in the development and validation of such approaches [12,13]. Moreover, many studies—particularly those with predictive objectives related to water resources under CC - rely on models in which hydrodynamic properties must be accurately defined [14,15]. With advances in computational material and methods, subsurface and groundwater flows are now widely represented in hydrological models using the Richards’ equation, particularly when applied at the catchment scale [14,15,16,17]. The hydrodynamic parameters used to describe the water retention and hydraulic conductivity curves of porous media can be estimated through pedotransfer functions, which are generally derived from particle-size analyses of soil samples [14,18,19]. Other approaches based on geophysical measurements can also be considered to characterize subsurface structures and properties in a less intrusive manner and/or to cover a larger spatial extent [20,21]. When soil moisture or drainage data are recorded in situ (e.g., TDR sensors and/or lysimeters), it becomes possible to relate hydrodynamic parameters to time series through an inverse modeling approach. Although such methods require numerical expertise and measurement efforts that can be complex due to field constraints, they nevertheless yield valuable and insightful results [19,22]. The presence of rock fragments in the soil is often neglected, as most hydrodynamic models are not designed to represent their influence. Nevertheless, introducing appropriate corrections may be essential to accurately evaluate water dynamics and storage in porous media [23,24].

A study conducted at the Strengbach catchment Observatory (Vosges Mountains, France) investigated the impact of monitoring frequency on the estimation of soil hydrodynamic parameters and proposed a consistent parameter set derived from matric potential measurements along a vertical soil profile [22]. Following severe forest dieback, largely attributed to pest outbreaks, the spruce stand at this site was completely clear-cut, and new monitoring systems were installed in 2022 on two more resilient plots.

The primary objective of the present study is to propose sets of hydrodynamic parameters suitable for hydrological modeling within the forested Strengbach catchment. These parameters are derived for two long-term monitoring plots, respectively dominated by spruce (Picea abies) and beech (Fagus sylvatica). The original results presented here are based on soil water content measurements obtained from two vertical profiles where soil samples were collected and analyzed. A dedicated methodology is described to refine the conventional approach based on PTFs by explicitly accounting for the influence of soil stoniness, thereby improving the representativeness of the derived hydrodynamic properties.

2. Materials and Methods

The following sections successively present the study site and the characteritics of the two plots studied. The instruments used for data acquisition are then described, along with the models employed to derive the hydrodynamic parameters. In sequence, the calculation of potential evapotranspiration and its implementation in the BILHYDAY model contribute to the estimation of actual evapotranspiration (AET). This variable is involved in defining the upper boundary condition used in the inverse modeling approach based on Richards’ equation, which is also presented. Finally, a specific section details the methodological aspects necessary to understand the steps involved in parameter estimation.

2.1. General Description of the Site

The Strengbach research catchment is located in the high-altitude commune of Aubure (Vosges Mountains, Haut-Rhin, France, see Figure 1). This experimental site covers an area of 0.8 km² and extends from 850 to 1150 m a.s.l. Approximately 80% of the catchment is forested, with about 76% consisting of a monospecific Norway spruce (Picea abies) plantation and the remaining area dominated by beech (Fagus sylvatica) [25]. From a pedological perspective, the catchment is characterized by acidic brown and ochreous podzolic soils developed on a calcium-poor granitic substrate [25].

Meteorological, hydrological, and geochemical variables of the Strengbach catchment have been continuously monitored since 1986 by the Hydro-Geochemical Environmental Observatory (OHGE; https://ohge.unistra.fr). The initial establishment of the monitoring network was motivated by research on acid rain and its effects on forest ecosystems. Current investigations now address a wider range of environmental issues, including forest decline associated with nutrient-poor soils [26,27], the impacts of climate change [28], and water resource management [29,30]. In addition, methodological developments are being pursued to improve the characterization of subsurface aquifers throughout multidisciplinary research [20,25,31].

2.2. Monitoring Devices and Associated Measurements

The Strengbach catchment is equipped with a meteorological station located in an open area at the top of the mountain (cf. Figure 1). Precipitation is measured using a tipping-bucket rain gauge (Précis Mécanique, model R5-302 b). Air temperature and relative humidity are recorded at a height of 1.5 m with a Vaisala HMP45C probe. Wind speed and direction are measured with an R.M. Young anemometer (model 05103), and incoming global radiation is monitored using a Kipp & Zonen pyranometer (model CM6B). All measurements are recorded by a Campbell CR1000X datalogger at a 10-minute time step.

For soil moisture monitoring, in the spruce plot (denoted JP in Figure 1), a soil pit was excavated to install Campbell TDR probes (model CS616) at depths of 13, 31, 52, 72, and 92 cm. In addition, four Campbell temperature probes (model 107) were installed at the first four depths to enable temperature correction. The sensor raw period (τ_raw, in μs) is corrected using a quadratic function of soil temperature (T_soil, in °C) defined in equation 1:

$${\tau }_{cor}={\tau }_{raw}+\left(20-T_{soil}\right)\times \left(5.26\times {10}^{-1}-5.2\times {10}^{-2}\times {\tau }_{raw}+1.36\times {10}^{-3}\times {\tau }_{raw}^2\right)$$

(1)

Following several laboratory tests conducted on soil samples, the standard quadratic calibration equation provided by the manufacturer was retained to derive the volumetric water content (θ, in m³/m³), as shown in equation 2:

$$\theta =-6.63\times {10}^{-2}-6.3\times {10}^{-3}\times \tau +7\times {10}^{-4}\times {\tau }^2$$

(2)

In the beech plot (denoted HET in Figure 1), a 50-cm Campbell SoilVue10 probe was installed, providing measurements of volumetric water content at depths of 5, 10, 20, 30, 40, and 50 cm. This sensor internally applies the correction and calibration equations and directly reports the volumetric water content. For both plots, WC data were recorded every 10 minutes on Campbell CR1000 (JP) and CR1000X (HET) dataloggers.

Notice that the SoilVue10 sensor does not require excavation of a pit but only the drilling of a hole using an auger with a slightly smaller diameter than that of the probe, into which it is then inserted. Owing to the pedological characteristics of the stand, a certain period was necessary for proper soil–probe contact to be established. A preliminary data processing step was applied to filter out measurements affected by preferential flow along the vertical probe and by water accumulation at the bottom of the hole. It should also be noted that this device operates at shallower depths and does not rely on exactly the same type of sensors as the installation on JP stand.

Although these data are primarily used as input for the inverse modeling process, the precipitation and soil moisture time series will be presented in the Results section to enhance the comprehension and assessment of the modeling outcomes. A complementary hydrological analysis will also be provided.

2.3. Vertical Profiles and Soil Properties

Undisturbed samples were collected in April 2022 from the JP plot at eleven depths (Figure 2a). Soil composition was determined using a laser diffraction particle sizing technique. The soil was mainly composed of particles larger than 2 μm, corresponding to sand and loam fractions, with the latter generally more abundant, averaging 53% compared to 35% for the sand fraction. For the HET plot (Figure 2b), data from [26] were used. They are based on the same analytical technique and show that the sand fraction is dominant, exceeding 59% across the depths analyzed, with a mean value of 65%, compared to 24% and 11% for the loam and clay fractions, respectively.

Sensor locations are shown in Figure 2a,b (red stars on the right side). Due to restrictions in operating the devices and the complexity of parameter estimation under natural conditions, the domain discretization was simplified. A step in our methodology described in section 2.5 will consist of selecting relatively homogeneous layers, each containing at least one sensor, and estimating their hydrodynamic parameters. Figure 2c,d show the pit excavated for sensor installation and sampling, with visible coarse elements that were not captured by the soil sampling cylinders of approximately 4 cm radius. The site was refilled while maintaining the original sequence of soil layers as much as possible, in order to monitor moisture dynamics in an environment minimally disturbed given the chosen methodology.

2.4. Numerical Methods

2.4.1. Potential and Actual Evapotranspiration

The choice of the potential evapotranspiration (PET) estimation method is crucial for accurately quantifying water fluxes in hydrological models [32]. To ensure consistency with previous studies conducted in the Strengbach watershed [22,33], the formulation derived from the Penman equation by Brochet and Gerbier [34] was adopted, as presented in Equation (3). Earlier studies have demonstrated that this approach accounts for the specific characteristics of the experimental site, particularly in the computation of net radiation in the radiative term (first component on the r.-h. s. term of Equation (3)) and in the formulation of the advective term (second component on r.-h. s. term of Equation (3)).

$${PET}_{BG}={\displaystyle \frac{1}{59}}{\displaystyle \frac{\Delta R_n}{\left(\Delta +\gamma \right)}}+{\displaystyle \frac{\gamma }{\left(\Delta +\gamma \right)}}0.26\times \left(1+0.4\times u_{10}\right){\delta }_e$$

(3)

PET_BG is given in mm.d^-1, Δ is the derivative of saturated vapor pressure versus temperature (mbar.K^-1) expressed as the air temperature Ta (°C) (see Equation (4)), R_n is the net radiation at the canopy surface (cal.cm^-2.d^-1), δ_e is the saturation vapor pressure deficit (mbar) expressed at the air temperature T_a (°C) and the relative humidity H_r (%) see Equation (5)), γ is the psychrometric constant (0.66 mbar.K^-1) and u₁₀ is the adjusted wind speed at 10 m height (m.s^-1).

$$\Delta ={\displaystyle \frac{25039}{{\left(237.3+T_a\right)}^2}}{.10}^{{\displaystyle \frac{7.5\times T_a}{237.3+T_a}}}$$

(4)

$${\delta }_e=\left(1-H_r\right)\times 6.107\times {10}^{{\displaystyle \frac{7.5\times T_a}{237.3+T_a}}}$$

(5)

Originally, the estimation of net radiation was derived from measurements of global radiation, complemented by a sunshine duration table used for potential corrections. Consequently, the preceding equations (Equations (3)–(5)) enabled the calculation of the daily PET_BG. Nevertheless, since meteorological measurements are available at higher frequency, PET can also be estimated at an hourly resolution (notice that daily and hourly values have been used in this study).

If PET represents the maximum evapotranspiration under well-watered conditions, the actual evapotranspiration (AET) accounts for water stress as well as precipitation intercepted and subsequently evaporated by vegetation foliage throughout the different seasons. An in-house code, BILHYDAY, based on a daily water mass balance and following a conceptual approach (see principles in [35]), was used to estimate AET according to Equation (6):

$$AET=Tr+In+SEvap$$

(6)

where Tr refers to transpiration by the canopy (mm.d^-1), In is the interception (mm.d^-1) and SEvap is the soil and understored evaporation (mm.d^-1).

Meteorological variables such as wind velocity, net radiation, air temperature, relative humidity, and recorded precipitation are provided to BILHYDAY as input data, along with site-specific parameters describing the canopy, plant interception, soil evaporation, and ecophysiological characteristics. PET_BG is then computed, and the three components of AET in Equation (6) are partly driven by the maximum leaf area index (LAI), a dimensionless variable that characterizes plant canopies and is defined as the projected leaf area per unit ground area. Since runoff has never been observed at the experimental site due to the continuous vegetation cover, the BILHYDAY model is well suited for estimating AET on a daily time step and has been previously calibrated on different plots of our research catchment [33,36]. The AET time series used to define the upper boundary condition will be described in the results section.

2.4.2. Vadose Zone Model

The adopted hydrological model is based on Richards’ equation (Equation 7) [37], which results from the combination of the mass conservation equation and the Darcy–Buckingham law, and describes water movement in a variably saturated porous medium.

$${\displaystyle \frac{\partial \theta }{\partial t}}+S{\displaystyle \frac{\theta }{{\theta }_s}}{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[K\left(h\right)\left({\displaystyle \frac{\partial h}{\partial z}}-1\right)\right]$$

(7)

where S is the storage coefficient (m^-1), θ is the volumetric water content (m³.m^-3), θ_s is the saturated water content (m³.m^-3), h is the pressure head (m), z is the vertical distance positive downward (m), t is the time (d) and K is the hydraulic conductivity (m.d^-1).

In addition, the standard Mualem–van Genuchten (MvG) model (Equations 8 and 9) [38,39] defines the relationships between negative pressure head, water content, and relative hydraulic conductivity (K_r, m.d^-1).

$$S_e={\displaystyle \frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}}={\left[1+{\left(\alpha \left|h\right|\right)}^n\right]}^{-m}$$

(8)

$$K\left(S_e\right)=K_{sat}\times K_r\left(S_e\right)=K_{sat}\times S_e^L{\left[1-{\left(1-{\left(S_e\right)}^{1/m}\right)}^m\right]}^2$$

(9)

where S_e is the effective saturation (-), θ_r is the residual water content (m³.m^-3), α is a parameter related to the mean pore size (m^-1), n is a parameter reflecting the uniformity of the pore size distribution (-), m = 1 – 1 / n is a parameter (-) [39], K_sat is the saturated hydraulic conductivity (m.d^-1) and L is a parameter related to the tortuosity chosen here equal to 0.5 (-) [38]. Notice that for saturated conditions (i.e. h ≥ 0), K = K_sat and S_e = 1.

The vertical profiles of forest soils contain a non-negligible proportion of stones and coarse fragments, which are generally not represented in undisturbed samples collected with cylinders and analyzed in the laboratory. To account for their effect on the evaluation of flow in the vadose zone, a corrective model for water content and hydraulic conductivity was implemented [23]. This model consists of two equations (10) and (11) for the water content and hydraulic conductivity respectively [23,24,40,41].

$${\theta }^b=\left(1-R_v\right){\theta }^f$$

(10)

$$K_{sat}^b=\left(1-a\times R_v\right)K_{sat}^f$$

(11)

In equations (10) and (11), the subscript b denotes the bulk value, whereas f refers to the fine soil fraction. The parameter a represents the hydraulic resistance of rock fragments to water flow and depends on their shape and size (-), while R_v denotes the relative volumetric fraction of stones (cm³.cm⁻³). The parameter a was set to 1.1 in accordance with [41] while the coefficient R_v is adjusted according to the plot and soil layer in question.

Hence, a set of 5 parameters (K_sat, θ_s, θ_r, n, α) is needed for each layer of the considered domain. Finally, the model requires specification of upper and lower boundary conditions throughout the simulation period.

2.4.3. Simulations Achieved

Simulations for the period 2022–2025 were performed using WaMoS-IPE-1D (Water Movement in Soil – Inverse Parameter Estimation – 1D) respectively for the spruce and beach plots. The numerical methods implemented in this Fortran code to solve Richards’ equation (Equation 7) include: (i) a finite element method for spatial discretization; (ii) a fully implicit backward Euler scheme for temporal integration; (iii) a variable switching approach employing modified Picard or Newton methods for linearization [42]; (iv) a heuristic time-stepping strategy based on iteration counts [43]; (v) a switching, flux-controlled procedure for the top boundary condition to maintain physically realistic results [44]; and (vi) a Levenberg–Marquardt iterative algorithm to solve the inverse problem [45]. A detailed description of the code is provided in Lehmann and Ackerer [46].

A uniform 1 cm mesh was adopted for the spatial discretization of both vertical profiles. Simulations were initialized with a specified hydrostatic pressure head profile derived from the lower water content measurement and allowed a sufficiently long warm-up period to ensure that initial conditions did not influence parameter estimation. An imposed flux boundary condition was applied at the soil surface, while the bottom boundary was defined as free drainage. As mentioned in section 2.4.1, the BILHYDAY code allows estimation of daily AET, which is then converted to hourly AET using equation (12). The top boundary condition is subsequently imposed at an hourly time step by subtracting the hourly AET to the rainfall.

$${AET}_h={\displaystyle \frac{PE{T}_{BG,h}}{PE{T}_{BG,day}}}{AET}_{day}$$

(12)

2.5. Methodology Investigated for the Model Calibration

The proposed methodology consists in the following steps:

• Definition of homogeneous soil layers

The first step consisted of subdividing the soil profile into homogeneous layers, each containing at least one sensor. Two different configurations to describe the 110 cm deep profile in the JP stand were set, while only one is used for the 100 cm deep HET profile.

2.

Initialization of hydrodynamic parameters

Field sampling and subsequent laboratory analyses enabled the determination of soil texture and composition for the different configurations described previously. Table 1 presents the configurations for the JP and HET stands. Following [22], initial estimates of the soil hydraulic parameters were obtained using the ROSETTA software based on PTFs [18], which derives parameters from measured soil texture data.

3.

Adjustments to account stoniness

The granulometric analyses performed on undisturbed soil samples were considered representative of the fine soil fraction. In contrast, both simulated and measured water contents represent the global bulk soil response. Consequently, the initial hydraulic parameters were corrected using the volumetric stone fraction (R_v) reported in Table 1 to account for the influence of stoniness on the soil water dynamics and the referring properties detailed in equations (10) and (11).

4.

Calibration and validation simulations

Given the available measurement ranges (which differ between the two sites), parameter inverse estimation was conducted using the WaMoS-IPE-1D model. The optimization procedure minimized an objective function defined as the mean squared error between observed and simulated water contents (θ_meas and θ_sim respectively):

$$Er={\displaystyle \sum _{n=1}^{N_{time}}{\displaystyle \sum _{i=1}^{N_{sens}}{\left({\theta }_{sim,i}^n-{\theta }_{meas,i}^n\right)}^2}}$$

(13)

where N_time is the number of observation times and N_sens denotes the number of layers instrumented with a sensor and used for the calibration.

Following [22], the inverse modeling approach focused on optimizing three key parameters: the saturated hydraulic conductivity (K_s), the shape parameter (α), and the pore-size distribution index (n). Hence, on the one hand, applying equation (10) is really important to get a proper estimation of θ_r and θ_s that are not estimated by our inverse model. On the other hand, K_s will be optimized and equation (11) is only useful for initialization. The quality of the calibrated parameter set was subsequently assessed using the Kling-Gupta Efficiency criterion [47] computed over the overall period of measurements according to equation (14):

$$KGE=1-\sqrt{{\left(r-1\right)}^2+{\left(\alpha -1\right)}^2+{\left(\beta -1\right)}^2}$$

(14)

where r refers to the correlation coefficient between simulated and observed WC, α to the ratio of their standard deviations, and β to the ratio of their mean values.

The Nash–Sutcliffe Efficiency (NSE) criterion [48] given in equation (13) can also be computed:

$$SE=1-{\displaystyle \frac{{\displaystyle \sum _{n=1}^{N_{time}}{\displaystyle \sum _{i=1}^{N_{sens}}{\left({\theta }_{sim,i}^n-{\theta }_{meas,i}^n\right)}^2}}}{{\displaystyle \sum _{n=1}^{N_{time}}{\displaystyle \sum _{i=1}^{N_{sens}}{\left(\overline{{\theta }_i}-{\theta }_{meas,i}^n\right)}^2}}}}$$

(15)

where $\overline{{\theta }_i}$ refers to the mean volumetric water content of the i^th layer over the period of concern. KGE and NSE values approaching 1 reflect an excellent agreement between simulated and observed data, while negative values indicate that the model performs poorer than simply using the mean of the observations as a predictor.

Depending on the plot considered (with different data availability), four calibration periods were selected:

• Calib.1: from 01 May 2022 to 30 September 2022 (no data available for HET),
• Calib.2: from 01 April 2023 to 31 December 2023,
• Calib.3: from 01 April 2024 to 31 December 2024,
• Calib.4: from 01 October 2024 to 30 September 2025.

One objective of the calibration step is to investigate the impact of the period used for calibration on the estimation of the parameters. The rest of the period is used for validation step.

3. Results

This section is dedicated to (i) the presentation of the meteorological time series and to the estimation of the mean site-scale potential evapotranspiration (PET), (ii) the presentation of the results obtained for each plot (JP and HET respectively) regarding the imposed upper boundary conditions and the outcomes of the parameter estimation using the inverse modeling approach.

3.1. Time Series of Climatic Forcing and PET

Meteorological data for the period of interest are presented in Figure 3. The hydrological year Y runs from October 1 of year Y to September 30 of year Y+1. Over the past 18 years, mean annual precipitation was 1223 mm, mean annual temperature 6.7°C, and total potential evapotranspiration (PET) 586 mm. The temporal evolution of daily precipitation and mean air temperature was used in the model described in section 2.4.1 to estimate potential evapotranspiration (PET). The four hydrological years considered exhibit substantial inter-annual variability relative to the 18-year climatological mean. The 2021 hydrological year, characterized by the summer drought of 2022, was the driest of the past 18 years (−20.7% rainfall, +18.8% PET relative to the long-term average), whereas 2023 was the wettest (+32% rainfall, -3.3% PET relative to the long-term average). Figures S1 and S2 illustrate this meteorological variability.

3.2. Results at JP Plot

Equation (11) was used to compute the AET at an hourly time step (based on BILHYDAY outputs). The resulting flux, used as a BC at the soil surface in the Wamos-IPE-1D model, corresponds to the difference between incoming precipitation and AET, as shown in Figure 4a. As commonly observed, the flux is positive, indicating a net meteoric water input to the soil. However, at the seasonal scale, the situation is more contrasted: AET accounts for between 44% and 86% (for the 2023 and 2024 hydrological years, respectively) of precipitation during the spring season (April–June), and between 57% and nearly 75% (for the 2023 and 2022 hydrological years, respectively) of precipitation during summer (July–September).

Accordingly, the temporal evolution of soil WC was analyzed at depths of 31, 52, and 72 cm, as illustrated in Figure 4b. This figure shows that the years 2022 and 2023 experienced marked and prolonged soil drying episodes, whereas in 2025, several rainfall events in August contributed to soil moisture recharge.

Equation (8) is applied to the parameters provided in Table 2, derived from ROSETTA, in order to obtain θ_r and θ_s values adjusted to the stoniness level of the soil profile. Then, the Levenberg–Marquardt algorithm implemented in WaMoS-IPE-1D is run to optimize the three parameters K_sat, α, and n. Table 3 summarizes the best set of parameters obtained for each layer, noting that our investigations led us to retain discretization JP2 and calibration period no. 3.

3.3. Results at HET Plot

Figure 5 a,b show the temporal evolution of the surface flux and of the soil WC at three depths (20, 30, and 50 cm) for the beech plot (HET). The descriptions are broadly similar to those provided for the spruce plot (JP). It is worth noting that sensor data could be reliably exploited starting from May 2023, although a few interruptions occurred thereafter. For this HET plot, the proportion of AET is approximately 10% lower than that observed for the spruce plot during spring, while the relative contribution during summer is slightly higher. The soil WC curves indicate that the spring and summer of 2023 (calendar year) correspond to the most pronounced soil drying period while the spring and summer of 2024 (calendar year) are the wettest recorded.

Following the same procedure as for the JP plot, the results obtained for the HET stand are presented in Table 4. The best set of parameters corresponds to calibration period no. 4.

4. Discussion

As highlighted in [22], the methodology focuses on optimizing three key parameters. To reduce the degrees of freedom of the inverse problem, the residual and saturated water contents are assigned fixed values. This approach is commonly applied in laboratory studies [49,50], whereas calibration using field data remains relatively scarce and requires additional supporting information [51].

Identifying parameter sets that are simultaneously optimal at both the global scale profile and for each individual layer is challenging, due to potential compensatory effects among flow processes in the different layers and, possibly, equifinality issues. Furthermore, the performance metrics are not always consistent across criteria. All results are provided in the supplementary information (Section 2 for the parameters and section 3 for the efficiency criterion). The optimal configuration retained for the JP plot (i.e., JP2-3) was obtained from calibration during period 3; except for layer 2, the efficiency indices are higher than those from the other configurations. For the HET plot, calibration period 4 (HET1-4) yielded the selected optimal parameter set. This configuration produced the best NSE values and offered the more balanced performance. In contrast, configuration 2 achieved the highest KGE values but the lowest NSE scores, leading us to discard it.

These findings highlight the need for caution when implementing an inversion technique. Calibration period 4 is the longest, with an average of 8,700 available measurements (across all three layers), whereas period 3 contains about 25% fewer data points. For the spruce plot, period 1 is the shortest and has the fewest soil moisture measurements. In the beech stand, period 1 contains no usable data, and period 2 is both the shortest and the least data-rich. These two calibration periods also exhibit the lowest average and minimum soil moisture values.

For the JP plot, the average soil moisture contents during calibration periods 3 and 4 are generally similar, while for the HET plot, they are slightly lower in 2025 than in 2024. In both plots, however, the 2024 hydrological year (calibration period 4) is characterized by lower minimum soil moisture values. The choice of calibration period therefore differs between the two plots. For the beech stand, the parameter sets HET1-4 and HET1-3 are relatively similar, whereas the set derived from calibration period 2 shows markedly different - and less efficient - parameter values. This discrepancy may be explained by particular environmental conditions or by reduced sensor reliability during that period, possibly due to imperfect contact between the sensor and the soil. The difficulty in accurately capturing the soil profile drying dynamics during the most intense episodes is likely due to the model applying transpiration only at the surface, whereas root water uptake should be represented over a root zone extending several tens of centimeters for both stands. In the supplementary information, Figure S7 shows the soil WC time series obtained with the HET1-2 model, which provides a better fit during the 2023 drying season but is no longer suitable for the rest of the simulation.

For the JP plot, the efficiency differences are more pronounced, making the choice of calibration period clearer. As with the beech stand, calibration period 2 results in a parameter set that is more distinct from the others, particularly regarding the properties of the successive soil layers. Finally, this calibration analysis supports the conclusions of [52], who emphasized the importance of considering both wet and dry conditions over a sufficiently long time span when estimating parameters based on soil moisture data. This study also reinforces the critical role of a subsequent validation step to ensure the robustness of the calibration results.

Furthermore, this study shows that direct simulations using parameter sets derived from PTFs (e.g., ROSETTA) lead to generally poorer performance, with KGE values of 0.74 for the JP plot and 0.26 for the HET plot, and NSE values of 0.15 and -2.75, respectively. This trend also holds when each soil layer is considered individually. The KGE result for the JP plot is somewhat misleading, as the agreement between simulated and observed soil moisture is not particularly good.

This work highlights that accounting for soil stoniness is essential to improve the model’s ability to reproduce field observations. Indeed, the parameters derived from laboratory measurements on soil samples differ substantially from the actual conditions encountered by the sensors in situ. Even when introducing an additional degree of freedom and including the estimation of the saturated water content (θ_s), the inversion results do not significantly improve model performance.

5. Conclusions

This study focused on the Strengbach catchment, which serves as the experimental site of the Hydro-Geochemical Environmental Observatory. Two forest plots, respectively composed of spruce (JP) and beech (HET), have been monitored for several years with soil water content sensors. However, the hydrodynamic parameters of the MvG model had not yet been determined for these stands. A methodology was therefore proposed to optimize, for each site, the set of parameters to be used in the Richards flow equation (Eq. 7). The results of this approach allowed us to identify optimal parameter sets (see Table 3 and Table 4, respectively) and to highlight several key points:

i. The PTFs method (here, ROSETTA was applied) proved functional but not optimal in this context.

ii. Taking into account soil stoniness is essential in this type of environment to improve flow model performance; the saturated water content is indeed strongly influenced by the presence of stones and gravel, confirming findings from earlier studies on the subject [40,41].

iii. Conventional performance criteria should be interpreted with caution, as they are not always consistent with one another. When calibration is performed over a limited period, it is essential to assess the robustness of the selected parameter set over an extended period.

iv. The selection of the calibration period is not straightforward; our study favored data representing a balance between dry and wet conditions (without focusing solely on drought periods) and covering a sufficiently long time span.

Our calibration targeted three parameters - α, K_s, and n - and it appeared that including the estimation of θ_s did not improve model performance, highlighting the complexity of characterizing the critical zone for hydrology modeling purposes. Several factors may explain the observed limitations and suggest possible improvements. On the one hand, other PTFs, potentially better suited to European soils (see [53]), or the integration of complementary information such as bulk density or field capacity, could enhance model efficiency. On the other hand, the model’s ability to capture low soil moisture values remains limited because transpiration is applied at the surface (i.e. throughout our BC), whereas it should be distributed over the rooting depth. Preferential flow paths, likely reinforced by the presence of coarse fragments, should also be considered. Finally, uncertainties remain in the estimation of environmental BV (precipitation and AET) as well as in sensor performance. The SoilVue10 probe, designed as a screw-type sensor, may experience partial disconnections from the surrounding soil, potentially affecting its readings; a longer probe length would probably be better suited. These aspects provide valuable matters for future research.