Three-way debate – measurements, pedotransfer functions, and inverse modeling – input strategy for soil hydraulic parameters in wheat simulations

1 Introduction

Yield variability in agricultural lands is expected to increase due to climate change, resulting in less reliable yields and jeopardizing food security (Asseng et al., 2013; Ray et al., 2015). The yield variability in cropland is driven by two dimensions: temporal variability, which reflects changes in weather conditions over time, and spatial variability, which refers to differences in soil and landscape characteristics across croplands (Basso et al., 2017). When analyzing yield variability in rainfed agriculture, it is crucial to consider in-season water that is temporarily present as rainfall and spatially stored in the soil. Influenced by pre-season and in-season rainfall, the amount of water retained in the soil profile is defined by its hydraulic properties, soil texture, structure, bulk density, porosity, and organic matter (Indoria et al., 2020). Plant-available water capacity refers to the amount of water that plants can access in the soil. The water retention capacity of the soil is defined by two key hydraulic parameters: field capacity and the permanent wilting point. Field capacity represents the maximum amount of water the soil can hold after excess water has drained away under gravity. On the other hand, the permanent wilting point indicates the minimum water content in the soil at which plant roots can no longer extract water. This point occurs when water is held so tightly in the soil that it becomes unavailable to plants. Accurately defining hydraulic parameters across diverse landscapes remains a significant challenge in crop production (Torres-Ruiz et al., 2024; Wesseling and Feddes, 2006). Pedotransfer functions (PTFs) are used as a cost-effective means to estimate soil hydraulic parameters from readily available data, such as soil texture. PTFs are developed using regression algorithms, where the correlation between soil properties and hydraulic parameters is determined and calibrated for a specific database (van Looy et al., 2017). The precision of empirical PTFs remains limited, even though PTFs have improved in recent decades; recalibration of PTFs is needed for new sites (Minasny and Hartemink, 2011; van Looy et al., 2017; Vereecken et al., 2016). Recent methodological advancements have integrated probabilistic methods to estimate soil water retention characteristics from limited measurements. Hierarchical Bayesian approaches have demonstrated significant utility in deriving site-specific retention curves from sparse datasets, enabling estimation of soil water retention behavior and quantification of associated uncertainty (Huang et al., 2025; Zhao et al., 2024). Although hierarchical Bayesian optimization is well established for inverse parameter estimation in broader hydrological and environmental disciplines, its explicit integration into process-based crop growth models for soil profile optimization remains scarce (Poudel et al., 2024).

By simulating the interactions among crops, management practices, and the environment, process-based crop models such as the Decision Support System for Agrotechnology Transfer (DSSAT) serve as analytical tools for evaluating rainfed agriculture (Araya et al., 2017). These models are used to support farm-level decision-making in marginal climates and to evaluate spatio-temporal cropping patterns. However, the predictive precision of these models is sensitive to the quality of the soil input data and the underlying water-balance logic (Rosa et al., 2012). Variability in soil properties affects water availability for plants, influencing the dynamic growth processes modeled in agroecosystem simulations (Paz et al., 1999). Because variations in plant-available water are widely recognized as a primary driver of yield gaps, inaccurate soil inputs in crop models propagate significant uncertainty into final yield predictions (Braga and Jones, 2004; Hoffmann et al., 2016; Wu et al., 2019). To capture these complex soil-plant-water interactions, process-based models rely on specialized sub-models. In the DSSAT-CSM framework, modules such as CERES-Wheat calculate soil water dynamics and root uptake using the one-dimensional Ritchie tipping-bucket water balance routine (Ritchie, 1998). The working mechanism of the tipping-bucket routine relies on field capacity and the permanent wilting point to define the plant-available water capacity. Within this routine, field capacity serves as the threshold regulating downward water movement through soil layers, while the permanent wilting point dictates the minimum threshold for root water extraction. If hydraulic parameters are not measured, DSSAT relies on internal PTFs to estimate them. Specifically, in the soil module, the PTF of Saxton et al. (1986) is used when soil texture data are provided, whereas the PTF of Rawls and Brakensiek (1985) is used when soil texture data and organic carbon are included as input variables. Model errors are diverse and numerous, reflecting the complexity of the soil-crop system, which is characterized by both known and unknown uncertainties (Boote et al., 2013). This is particularly evident in estimating soil hydraulic parameters, where researchers often face a ‘swamp of complete confusion’ due to high variability and a lack of consensus among existing estimation methods (Gijsman et al., 2002). Within the CERES-Wheat model, this confusion directly increases uncertainty when the in-situ field capacity is simply equated to laboratory-measured field capacity (assessed at a matric potential of -33 kPa). While modern laboratory protocols appropriately rely on undisturbed soil cores to preserve local structure, these isolated samples inherently fail to account for the natural layering of soil horizons, broader structural dynamics, and profile-scale drainage restrictions present in the field, which can lead to significant simulation errors (Dane and Hopmans, 2002). The lower threshold of plant-extractable water in the field is often treated as synonymous with the laboratory-measured permanent wilting point (-1500 kPa). To navigate these uncertainties and reconcile theoretical tipping-bucket logic with field observations, systematic optimization is essential (van Looy et al., 2017; Vereecken et al., 2016).

Inverse modeling has emerged as a robust approach for calibrating soil hydraulic parameters within process-based models (He et al., 2022; Hoffmann et al., 2016; Trenz et al., 2024). Parameter optimization is needed to reduce prediction error and variability in soil-crop models, but the effectiveness of such inverse modeling depends on the optimization strategy to estimate new parameters (Wallach et al., 2024). The techniques range from straightforward systematic grid search to more complex heuristic approaches such as genetic algorithms or simulated annealing (Press et al., 1992). The selection of a specific method involves an inherent trade-off among computational efficiency, problem complexity, and the guarantee of reaching a global rather than a local optimum (Mayer et al., 2001). However, inverse parameter optimization carries a significant risk of equifinality, in which multiple parameter combinations yield an excellent statistical fit yet fail to represent a realistic physical soil-water environment. This issue primarily arises from the fundamental structure of crop models and pedotransfer functions (PTFs). Within frameworks like CERES-Wheat, certain model inputs exhibit inter-dependence or function as empirical coefficients without a clear, measurable physical meaning. For instance, the whole-profile drainage rate coefficient and layer-specific root growth factors serve as conceptual tuning variables rather than directly measurable field properties (He et al., 2022; Jones et al., 2003). Furthermore, hydraulic parameters such as field capacity and permanent wilting points are inherently interdependent in defining plant-available water capacity; consequently, multiple contrasting pairs can mathematically yield the same simulated root water uptake and crop yield, leading to equifinality (Beven and Freer, 2001). Relying on generalized PTFs can occasionally exacerbate this issue, as their underlying empirical relationships may struggle to separate correlated soil traits. To mitigate equifinality and ensure the physical validity of the simulated environment, a critical evaluation of both the model structure and the applied PTFs is necessary. Consequently, parameter optimization must be approached through a systematic, stepwise calibration strategy. By first anchoring and adjusting the explicitly measurable, physically based parameters before fine-tuning the conceptual or non-physical variables, the optimization space is effectively constrained, facilitating the derivation of realistic parameter values. To ensure the models simulate the right results for the right reasons, the objective function must be conditioned by the appropriate observed variables (Beven and Binley, 1992).

The objective of this study is to systematically reduce simulation uncertainty in the DSSAT-CERES-Wheat model by evaluating the stepwise transition from directly measured to empirically estimated, and ultimately to inversely optimized soil hydraulic limits. This evaluation is grounded in a comprehensive 3-yr dataset featuring high-resolution in situ measurements of volumetric soil water content and matric potential. The objectives are: (a) Validate the predictive accuracy of Saxton and Rawls default PTFs by comparing predicted retention curves against lab-measured and sensor-derived data. (b) Analyze the model’s response to variations in field capacity and permanent wilting point by simulating soil water content across distinct soil horizons (15 and 50 cm). (c) Benchmark the efficiency of Bayesian versus Grid search optimization algorithms within the SPO software plug-in to determine the most robust approach for inverse modeling if wheat growth variables can effectively calibrate soil field capacity and permanent wilting point parameters. (d) Appraise the integrated soil-water balance to confirm that the optimization algorithms and their predictions remain within physically realistic hydrological dynamics.

2 Materials and methods

2.1 Description of the study site

A 3-year on-farm field study (2022, 2023, 2025) with wheat (Triticum spelta L. cv. Bauländer Spelz) was carried out at the field ‘In der Klinge’, located in Boxberg, northern Baden-Württemberg (Germany), at 358 m above sea level (49°45N, 9°65E), which was part of the research project CHARGE (Field investigation on the impact of high voltage direct current (HVDC) underground cables on soils and crop production, PTKA/KIT L7521101); the experimental data of the 2022 and 2023 seasons were published in Trenz et al. (2025) and further information can be found there. The year 2024 was excluded from the modeling dataset because maize was cultivated, thereby maintaining experimental consistency by evaluating the same wheat cultivar across all analyzed seasons. The climatic conditions for the past 25 years (2001-2025) are characterized by a mean annual precipitation of 751 mm, a mean air temperature of 9.8 °C, and global radiation of 978 kWh m^-2. Over the experimental years, the mean annual temperature exceeded the 25-year average (11.1 °C in 2022, 11 °C in 2023, and 10.2 °C in 2025). According to the IUSS Working Group WRB (2022), the soil profile at the experimental site is classified as a vertical Cambisol with a dominant silt loam texture. Table 1 details the soil profile, illustrating its five natural soil horizons (Ap1, Ap2, Bw, Big1, and Big2), alongside the soil texture (clay, sand, and silt), and organic carbon content (C_org). Soil texture was determined using the pipette and sedimentation method in accordance with ISO 11277:2020, 2002 and the KA5 German classification fractions (Sand: 2.0–0.063 mm; Silt: 0.063–0.002 mm; Clay:<0.002 mm) (Boden, 2005). To account for differences in particle-size limits between classification systems, these data were transformed to the USDA scale following the methodology of Rawls et al. (1982). This conversion was performed via a log-linear transformation using ‘The Soil Texture Wizard’ R package (Moeys, 2010).

Table 1. Soil profile characterization for the five natural horizons at the Boxberg site, including soil texture (sand, silt, and clay) after KA 5 (Boden, 2005) and USDA after Rawls et al. (1982), C_org (%), and USDA texture classification, is shown.

Across all three experimental seasons, wheat was drill-seeded using an Accord pneumatic airflow system at a depth of 5 cm, with an inter-row distance of 12.5 cm and a seeding rate of 150 kg ha^-1 (290 seeds m^-2), on November 6, 2021; November 12, 2022; and November 9, 2024. On March 3, 2022, 20 kg N ha^-1 was fertilized as calcium ammonium nitrate (CAN, 27% N). Wheat was harvested on 18^th July 2022. On April 4, 2023, 50 kg N ha^-1 was fertilized as CAN, and wheat was harvested on the 18^th of July, 2023. In 2025, CAN was fertilized on the 9^th of March 2025 with 45 kg N ha^-1, and on the 11^th of April 2025 with 25 kg N ha^-1. Wheat was harvested on the 14^th of July.

2.2 Data collection and analysis

2.2.1 Plant data

The experimental field covers 0.8 ha and was divided into four columns, each consisting of six aligned, site-specific plots, for a total of 24 plots (Figure 1). Each site-specific plot measured 9 × 9 m. Plant sampling involved destructive biomass cuts conducted five times throughout the growing season. Biomass cuts were randomly taken with a size of 0.0625 m² (0.25 m × 0.25 m) within each plot. Wheat development ratings were conducted on the Zadoks scale (Zadok et al., 1974). The first three destructive cuts occurred during the vegetative growth phase of wheat (BBCH 31, BBCH 35) and at the flowering stage (BBCH 65). Above-ground biomass was collected and separated into plant organs, including leaves, stems, and ears. The fourth biomass cut was performed during the generative phase, between grain filling and ripening (BBCH 75). The final biomass cut was conducted at maturity. After organ separation, the samples were dried at 60 °C to a constant weight to determine dry matter weight. Leaf Area Index (LAI) was measured non-destructively within the plots using an LAI-2200C device (LI-COR, Lincoln, USA).

2.2.2 Soil data

The soil properties C_org, pH, and the sand, silt, and clay fractions were determined from bag samples collected from two profile pits. For all subsequent procedures, the samples were first air-dried at 40 °C and then sieved to 2 mm. CO₃^- and C_org were determined using a SoliTOC Cube (Elementar Analysensysteme GmbH, Germany). The pH was measured potentiometrically in a 0.01 M CaCl₂ solution and a Mettler Toledo meter (SevenEasy™ S20) equipped with an In-Lab Expert Pro electrode.

To determine the soil hydraulic characteristics, five standard pF-θ data points were measured representing saturation and matric potentials of pF 1.0, 1.8, 2.5, and 4.2 (corresponding to -1 kPa, -6.3 kPa, -33 kPa, and -1500 kPa, respectively). The pF value is defined as the base-10 logarithm of the absolute soil matric potential, expressed in centimeters of water, thereby rendering the pF index dimensionless. The corresponding volumetric soil water content θ is expressed in units of cm³ cm^-3. Undisturbed soil cylinders (100 cm³ volume; inner diameter: 5.64 cm, height: 4 cm) were placed in a vessel lined with a shallow layer of coarse sand. Over a three-day period, a peristaltic pump was used to gradually saturate the samples from the base upward. The saturation solution consisted of 5 mM CaSO₄ and 2 mM thymol to maintain chemical stability and inhibit microbial growth. Following saturation, specific matric potentials were established: pF levels 1.0 and 1.8 were attained using a hanging water column (tension method), while pF levels 2.5 and 4.2 were achieved using the overpressure method (Schlichting et al., 1995). The pF 4.2 measurement was conducted on disturbed soil suspensions in 4.5 cm by 1 cm PVC rings. To determine the volumetric soil water content (θ) at each pF step, the samples were weighed upon reaching hydraulic equilibrium. After the final pressure step, the samples were oven-dried to a constant weight to determine the total dry soil mass. The volumetric water content (θ) for each pressure step was subsequently calculated based on the difference between the wet and dry mass relative to the 100 cm³ sample volume. Bulk density was simultaneously computed from this final dry soil mass. Crucially, these standard laboratory measurements provided baseline estimates for the crop model’s soil-hydraulic boundaries. Aligning with standard DSSAT conventions, the laboratory water content at pF 2.5 served as the initial proxy for field capacity, while the measured water content at pF 4.2 served as the baseline estimate for the permanent wilting point. The retention curves were fitted, and estimated data pairs of matric potential (pF) and volumetric water content (cm³ cm^-3) were obtained using the van Genuchten (1980) model described in Equation 1:

(1) θ(h)=θr+θs−θr(1+(α|h|)n)m

θ(h) is the volumetric water content at a given pressure head (h). The parameters θ_s and θ_r represent the saturated and residual volumetric water contents, respectively. The empirical parameters α, n, and m determine the shape of the retention curve. Specifically, α acts as a scaling parameter relative to the matric potential axis, while n is a measure of the pore-size distribution (van Lier and Pinheiro, 2018). In this study, the Mualem constraint m = 1-1/n was assumed to reduce the number of independent parameters and allow for the prediction of unsaturated hydraulic conductivity.

Soil nitrogen levels (N_min) were measured at three depths (0–30, 30–60, and 60–90 cm) before sowing and after harvest in two adjacent plots, for a total of 12 replicates. The N_min content was determined following the VDLUFA (1991) method: homogenized soil samples were extracted with a 0.01 mol L^-1 CaCl₂ solution at a 1:10 soil-to-solution ratio and shaken for 1.5 h. The resulting extracts were analyzed to quantify nitrate-nitrogen and ammonium-nitrogen across all soil depths. Volumetric soil water content and matric potential were monitored at depths of 15, 50, and 130 cm using TEROS 11 (± 3 vol.% accuracy) and TEROS 21 (± 10% of the reading + 20 hPa accuracy) sensors, respectively (Meter Group GmbH, Munich, Germany). The TEROS 11 sensors were distributed across six plots, while the TEROS 21 sensors were installed in three plots (Figure 1). The TEROS 21 operates as a dielectric water potential sensor, measuring the dielectric permittivity of a calibrated, porous ceramic disc in hydraulic equilibrium with the surrounding soil, and reports values on the dimensionless pF scale. Regarding the installation timeline, subsoil sensors (50 and 130 cm) were permanently installed and provided continuous data for 2022 and 2023; however, a data transmission failure precluded subsoil recordings in 2025. To comply with farm soil conservation practices, topsoil sensors (15 cm) were installed each year after spring field operations temporarily, capturing data throughout the vegetative periods of all three seasons.

2.2.3 Data analysis

All in-season soil and plant measurements were analyzed using a linear mixed-effects model via the GLIMMIX procedure in SAS software (Version 9.4, SAS Institute Inc., Cary, NC, USA) to account for the split-plot design. The model considers block, column, and plot effects embedded within the blocks. The model was utilized to extract the estimated marginal means (Least Squares Means) for the undisturbed treatment across the 24 plots. The model for single-season analysis is described in Equation 2 as:

(2) yno=μ+sn+bno+eno

Where y_no is the observation in column (or main-plot) n and block o, and μ is the intercept. The terms s_n and b_no are the random main-plot and block effects of column n and block o, respectively. The term e_no denotes the residual error, which accounts for the variance of individual plots within blocks. These estimated means represent the field-specific baseline, which was subsequently used as input for the modeling study. Deriving these field-specific means via the mixed model provides considerably more accurate estimates than simple arithmetic averaging, as it properly accounts for the variance introduced by the spatial and block structures (Littell et al., 2006).

2.3 The DSSAT crop simulation model and soil water balance

Process-based agroecosystem models are fundamental mathematical tools that simulate the complex soil-plant-atmosphere continuum, calculating daily mass and energy balances to predict crop growth and resource use. The Decision Support System for Agrotechnology Transfer (DSSAT v4.8.5) is a comprehensive software application that integrates crop simulation models for over 40 crops, thereby enhancing predictive capabilities (Boote et al., 2025; Hoogenboom et al., 2019). The DSSAT CSM-CERES-Wheat framework was specifically selected for this study. Because the overarching objective of this work is to reduce uncertainty in yield predictions due to inaccurate soil inputs, CERES-Wheat provides an ideal testing environment. Its soil-water balance is sensitive to defined hydraulic boundaries, making it well-suited for evaluating how optimized parameterization reduces model uncertainty. The model predicts wheat production by accounting for water and nitrogen balances using cultivar-specific coefficients. The CSM-CERES framework estimates potential dry matter production using a radiation-use-efficiency approach and adjusts the results based on soil state variables such as water and nitrogen availability. Furthermore, the accumulation of degree days, photoperiod sensitivity, and vernalization requirements fundamentally govern the simulated crop development stages.

Soil water dynamics within the DSSAT-CSM-CERES framework are controlled by a one-dimensional, layer-wise tipping-bucket water balance module adapted from Ritchie (1972, 1981). To account for discontinuities in water movement across the varying soil horizons, the model uses a multilayered soil profile structure. The simulated profile is not treated as a single homogeneous unit; rather, it is discretized into separate computational layers. Each layer is parameterized with its own specific retention limits and conductivity parameters.

Consequently, the model’s water balance module calculates vertical water fluxes based on the distinct hydraulic properties of each individual layer, capturing the natural discontinuities in water movement across the entire profile. Functional soil water storage in each discrete layer is bounded by three critical hydraulic parameters: the permanent wilting point (referring to -1500 kPa), field capacity (referring to -33 kPa), and the saturated water content (bounded by total soil porosity). In this study, standard soil physics terminology is used; the terms field capacity and permanent wilting point correspond directly to the DSSAT model inputs for the drained upper limit (volumetric water content at -33 kPa) and the lower limit (volumetric water content at -1500 kPa), respectively. While these hydraulic thresholds can be measured directly in the laboratory, such data are not always available. Therefore, in the absence of direct hydraulic measurements, these parameters are estimated from a minimal set of basic, measured soil data using established pedotransfer functions (PTFs). Depending on the input data, the DSSAT Soil Profile Creator (SBuild) uses different PTFs to calculate soil hydraulic parameters at saturation, at field capacity, and at the permanent wilting point. In SBuild, it is selected between two PTFs based on the availability of soil organic carbon (SOC =C_org) data: when C_org is present for all horizons, the Rawls et al. (1982) PTF is used; when C_org is missing in any horizon, the Saxton et al. (1986) PTF is used as a fallback. Both PTFs estimate soil water content at saturation (SAT) as 95% of total porosity (POR) (SAT = 0.95 × POR), where porosity is derived from an organically adjusted particle density, following Baumer and Rice (1988). The PTFs diverge in their estimates of the field capacity and permanent wilting point. Rawls and Brakensiek (1985) use empirical regressions at 1500 kPa and 33 kPa, respectively, incorporating clay, sand, and soil organic matter (SOM) as predictors with SOC converted to SOM using a factor; Saxton et al. (1986) derive field capacity and permanent wilting points from texture alone, substituting a fixed SOM value of 0.66% when organic matter data are not available. The Rawls PTF yields systematically higher estimates of field capacity- and permanent wilting point values in organically enriched topsoil horizons, with differences diminishing with depth as SOC decreases. Both PTFs calculate based on layer-wise soil texture (sand and clay mass fractions), organic matter content (OM), and soil bulk density ρ_b the hydraulic limits—permanent wilting point (θ₁₅₀₀), field capacity (θ₃₃), and saturation (θ_s) (Equations 3–5):

(3) θ1500=f(Sand,Clay,OM)

(4) θ33=f(Sand,Clay,OM,ρb)

(5) SAT=f(θ33,Sand,ρb)

Where θ₁₅₀₀ and θ₃₃ represent volumetric soil water contents at -1500 kPa and -33 kPa. Infiltration is calculated by subtracting surface runoff, estimated using the USDA Soil Conservation Service Curve Number (CN) method, from the sum of total daily rainfall and irrigation (Soil Conservation Service, 1972).

Reference evapotranspiration (ET₀) is estimated using the Priestley–Taylor method and partitioned into maximum plant transpiration (E_p) and maximum soil evaporation (E_s) based on LAI. Actual soil evaporation is simulated using the classical Ritchie-CERES model (Ray et al., 2015; Ritchie, 1972). In the first stage, evaporation occurs when the surface soil is wet, proceeding at the maximum energy-limited rate. Once cumulative evaporation exceeds a soil-specific threshold, the second stage begins, in which actual soil evaporation becomes soil-limited and declines as the square root of time since drying started. Downward drainage occurs sequentially between layers whenever volumetric water content exceeds field capacity, governed by a layer-specific drainage fraction (SLDR, day^-¹).

2.3.1 The soil profile optimizer: grid search and Bayesian optimization options

Trenz et al. (2024) introduced a generic optimization approach for soil profile parameters, implemented using an external software plug-in called Soil Profile Optimizer (SPO), available on GitHub (https://github.com/memicemir/SPO) and compatible with DSSAT v4.8.5. The SPO is a unifying inverse modeling approach that optimizes the objective function using parameter search techniques based on available crop measurements, such as LAI, grain weight, and above-ground biomass, as well as soil measurements, such as layer-based soil water content. The SPO is based on an inverse modeling approach that employs a combination of Grid search and Bayesian optimization techniques to iteratively calibrate specific soil hydraulic parameters from observed data. The SPO adapts selected soil profile parameters within a defined range by minimizing the difference between simulated and observed target variables, as measured by the nRMSE. Any in-season data observed can serve as target variables in the software, including yield, above-ground biomass, LAI, and soil water content (Trenz et al., 2024). The selection of specific optimization methods involves a trade-off between computational simplicity and search sophistication (Braga and Jones, 2004). These techniques range from straightforward Grid search, which systematically evaluates a predefined range of values, to more complex heuristic approaches, such as genetic algorithms designed to navigate high-dimensional parameter spaces more efficiently. Therefore, the current version of the SPO was refined, and the revised optimization procedure now allows the user to select between two optimization methods: Grid search and Bayesian algorithms. The Grid search method performs an exhaustive evaluation of the parameter space by discretizing the defined bounds of the soil hydraulic properties – the permanent wilting point and the field capacity across all soil layers – into a structured grid of fixed step intervals. Every possible parameter combination within the grid is evaluated by running DSSAT-CSM and computing the nRMSE between simulated and observed data, and the combination yielding the lowest error is selected as the optimal solution. While this approach is fully transparent and deterministic, and guarantees that no region of the defined search space is overlooked (Bergstra and Bengio, 2012), its computational complexity grows exponentially with the number of parameters and the grid resolution, making it less practical for high-dimensional calibration problems. The Bayesian optimization method is implemented using the open-source Python library developed by Nogueira (2014). The parameter search space was constrained to ±15% of the initial baseline soil values to ensure physical realism. The underlying surrogate model was a Gaussian Process Regressor using a Matérn kernel (ν = 2.5), which provides a flexible, smooth approximation of the DSSAT parameter response surface. The search process was guided by the Upper Confidence Bound acquisition function (κ=2.576) to balance exploration of uncertain regions of the parameter space with exploitation of regions with known low nRMSE. Each optimization run consisted of an initial random exploration phase of 10 points, followed by 100 Bayesian iterations. To assess stochastic robustness and account for the influence of the initial random sampling on convergence, the entire calibration pipeline was structured to perform repeated evaluations. The Bayesian optimization process was executed independently across three distinct random initialization seeds (42, 123, and 999). The final calibrated soil parameters reported for the model were drawn from the independent run that achieved the lowest global nRMSE, ensuring convergence stability and reproducibility.

2.3.2 DSSAT model setup, model calibration, and simulation settings

The collected soil and crop growth data from 2022, 2023, and 2025 were entered into the DSSAT v4.8.5 platform. Daily weather inputs, including minimum and maximum temperatures, solar radiation, and rainfall, were sourced from a weather station on the experimental site. These variables met the model’s climatic input requirements, as evapotranspiration was computed using the Ritchie method. Management practices, including planting and fertilization schedules, were recorded for three years and incorporated into the experimental files. For measured input data, the estimated means of a mixed model are used. The different simulation approaches in this study stem from varying soil input data, either arising from the initial model mechanism or generated during the optimization procedure. Table 2 provides an overview of the five simulation approaches, each based on a different parameterization of soil hydraulic parameters, applied in this study to evaluate their performance. DSSAT automatically divides the soil profile into standard depth increments, or can be configured to preserve the horizons exactly as measured in the field, rather than treating the profile as a single homogeneous unit or using arbitrary computational layers. Based on the in-field observations, the soil profile contained five natural soil horizons (Ap1, Ap2, Bw, Big1, and Big2) and was divided into nine soil layers. To maintain numerical stability in water balance calculations and accurately capture soil evaporation, DSSAT guidelines recommend topsoil layers of 5–10 cm, while deeper subsoil layers are limited to a maximum thickness of 30 cm (Ritchie, 1998). If a measured field horizon exceeds this 30 cm threshold, a subdivision into thinner, equally proportioned layers that retain the exact same physical and hydraulic properties as the original parent horizon is conducted. The simulation approach PTF (Rawls) was used in this study for initial model parametrization, calibration, and evaluation. The soil input files contain mean soil profile information derived from observed soil texture and C_org data. Hydraulic parameters (saturation, field capacity, and permanent wilting point), along with a comprehensive set of additional soil physical and hydraulic parameters, were estimated using the standard SBuild procedure. A similar procedure was followed for the PTF approach (Saxton); however, the soil input data included texture information but lacked C_org, which was added to SBuild. In the Lab-Measured approach, texture data, C_org, and measured soil water content for the field capacity (pF 2.5) and permanent wilting point (pF 4.2) were entered rather than calculated with SBuild. The optimization approaches were generally distinguished into Bayesian and Grid search. For both optimization approaches, the soil profile of the PTF (Saxton) approach was used as the starting point. In both approaches, the three target variables, soil water content (cm³ cm^-3) at 15 cm, grain weight (kg ha^-1), and above-ground biomass (kg ha^-1), were selected, and as a statistical measure, the nRMSE was used to minimize the difference between observed and simulated values. In both approaches, the search space buffer for the field capacity and permanent wilting points was set to ±15% for Bayesian and ±20% for Grid search. The cultivar coefficients were calibrated using the data from the 2023 season, while the 2022 and 2025 seasons were reserved for evaluation. The cultivar coefficients were calibrated manually, following the protocol established by Boote (1999). The calibration followed a sequential logic: first, wheat phenology was calibrated to anchor the growing season. Second, the vegetative phase—including leaf weight, stem weight, and LAI—was calibrated from emergence to anthesis to ensure the model accurately captured canopy architecture. Finally, the reproductive phase was refined using grain weight, tiller number, thousand-kernel weight, and grain number as key state variables. All simulations were conducted with the PTF (Rawls). In Supplementary Figure 1, the comparison of observed and simulated data for six growth variables (LAI, leaf weight, stem weight, above-ground biomass, kernel number, and grain weight), along with the d-statistics and the NSE, demonstrates the model calibration results for the PTF (Rawls) treatment. The calibration was conducted for the 2023 growing season, which benefited from more favorable conditions, including higher precipitation and higher N_min, resulting in higher yield potential. Although the NSE for leaf weight was low (-1.17), the overall calibration for the in-season parameters was excellent. Across all six growth variables, a d-statistic ≥ 0.76 was observed, indicating a high overall goodness of fit.

Table 2. Overview of the five simulation approaches based on different soil hydraulic parameterizations applied in this study to evaluate their model performance.

2.4 Statistical evaluation of simulation results

The statistical method used in this study to assess the agreement between observed and simulated data and to evaluate model performance is based on the following statistical measures. The index of agreement d statistic is a dimensionless sum-of-squares-based measure (0 ≤ d ≤1) of the extent of deviation from zero originally reported by Willmott (1982). O_i denotes the observed values, P_i denotes the simulated values, n denotes the number of observations, and Ō denotes the mean of the i^th observation. Therefore, d-statistics can be used for comparing multiple output variables (Equation 6):

(6) d=1−∑1=1n(Pi−Oi)2∑i=1n(|Pi−O¯|+|Oi−O¯|)2

The Nash-Sutcliffe model efficiency coefficient (NSE) was calculated alongside the index of agreement (Nash and Sutcliffe, 1970). The NSE determines the relative magnitude of the residual variance compared to the measured data variance, as shown in Equation 7:

(7) NSE=1−∑i=1n(Oi−Pi)2∑i=1n(Oi−O¯)2

where O_i is the i-th observed value, P_i is the i-th simulated value, Ō is the mean of the observed values, and n is the total number of observations. NSE values range from −∞ to 1.0, with an NSE of 1.0 indicating perfect agreement between simulated and observed values. An NSE of 0 indicates that the model predictions are as accurate as the mean of the observed data, while an NSE less than zero (NSE<0) indicates that the observed mean is a better predictor than the model.

The root mean square error (RMSE) as defined in Equation 8 in model calibration and validation for assessing the deviation (y-x) between simulated and observed values, providing a consistent measure of that deviation (Loague and Green, 1991).

(8) RMSE=[∑i=1n(Pi−Oi)2n]0.5

Additionally, nRMSE was used to evaluate the model. The nRMSE was calculated using Equation 9 and can serve as a relative measure for comparing state variables with different units (Loague and Green, 1991):

(9) nRMSE=RMSEO¯×100

Curve fitting was performed using the Levenberg-Marquardt nonlinear least squares algorithm ('nls.lm function'), which is robust to poor initial parameter estimates (Elzhov et al., 2022). Fitting was conducted separately for each soil layer and treatment. A two-attempt strategy was applied: the first attempt used fixed starting values and bounds based on typical ranges of soil hydraulic parameters. If the first attempt either failed to converge or produced physically implausible results (e.g., θs ≤ θr, or θs outside the observed data range), it was rejected, and a second attempt was made using data-driven starting values and bounds derived directly from the observed water content range of each layer. Goodness-of-fit and model performance across treatments were evaluated using the coefficient of determination (R²), calculated using Equation 10:

(10) R2=1−∑i=1N(θobs,i−θsim,i)2∑i=1N(θobs,i−θ¯obs)2

where θobs,i and θsim,i represent the measured and model-predicted volumetric water contents at i-th matric potential step, and θ¯obs is the mean observed volumetric water content.

3 Results

3.1 Seasonal soil matric potential

The rainfall amount and pattern during the vegetation period varied between 2022, 2023, and 2025. Considering the total precipitation from seeding to maturity (10^th of July), the 2022 season accumulated 452 mm, the 2023 season 508 mm, and the 2025 season 417 mm. Figure 2 shows the matric potential (pF) at 15 cm depth across the growing seasons (2022, 2023, and 2025) in combination with the daily rainfall. The red-shaded band marks a critical stress threshold zone (pF 2.7-3.0), above which plants normally experience water deficits. In 2022, the matric potential rose into the stress-threshold (pF > 2.7) zone between anthesis and maturity. In 2023, the pF values remained<2.0, with notable rainfall events. In 2025, a steady increase in pF values was observed from March to July, ending at pF 2.65 at maturity. As indicated by the error bars, measurement uncertainty was observed in the topsoil, likely reflecting rapid environmental dynamics and fluctuating soil-water contact. In Figure 3, soil matric potential (pF) at depths of 50 and 130 cm across the growing seasons of 2022 and 2023 is shown, with daily rainfall indicated by blue bars on the secondary y-axis. In 2022, at 50 cm, the matric potential exceeded the threshold band six days after anthesis (15^th of June), peaking above pF 3.0, and then dropped after harvest in October, with no significant rainfall during this period. In 2023, at a depth of 50 cm, matric potential increased to pF 2.8 on the 26^th of June, ten days after anthesis, but remained below the threshold band, aided by significant rainfall events between anthesis and maturity. At 130 cm, the highest matric potential in 2022 was pF 2.3, reached on the 7^th of July, while in 2023, a lower peak of pF 1.8 was observed on the 10^th of July, indicating reduced stress with soil depth. Comparing the two years, 2022 showed consistently higher matric potential values at both depths, reflecting higher overall stress. In both years, a more pronounced and faster increase in matric potential was observed at 50 cm than at 130 cm.

3.2 Measured vs calculated hydraulic parameters and their impact on soil water simulations

Figure 4 shows the fitted van Genuchten retention curves at three soil depths (15 cm, 50 cm, and 130 cm) for the four data sources: Sensor—in-situ field measurements (green), Lab-Measured—laboratory-derived retention curve (black), PTF(Saxton)—pedotransfer function estimate using the Saxton approach (orange), and PTF(Rawls)—pedotransfer function estimate using the Rawls model (yellow), plotting volumetric water content on the y-axis against matric potential on the x-axis. Goodness-of-fit was evaluated by calculating R² from non-linear regression models (van Genuchten fits) comparing volumetric water contents from lab measurements, pedotransfer functions (PTFs), and field sensors across three soil depths. Across all data sources, model performance was highest at 15 cm and 50 cm (R² = 0.89- 0.99), whereas the PTFs showed a drop in accuracy at 130 cm (R2 = 0.57- 0.72) compared to Lab-Measured (R²≥0.94). Additional information is given in Supplementary Table 1; soil water content at pF 2.5 (field capacity), pF 4.2 (permanent wilting point), and plant-available water capacity are listed for the Lab-measured, Sensor, PTF (Rawls), and PTF (Saxton), and the Percentage Deviation (PD%) for the plant-available water capacity was calculated. At 15 cm depth, sensor data cluster around a soil water content of 0.30-0.35 cm³ cm^-3, generally below the Lab-Measured curve, whereas the PTF estimates tend to overestimate water content compared to sensor readings, with the highest soil water content (0.517 cm³ cm^-3) at pF 0 for the PTF(Rawls). The Lab-measured water content at the field capacity was 0.303 cm³ cm^-3, and at the permanent wilting point was 0.125 cm³ cm^-3, resulting in a plant-available water capacity of 0.179 cm³ cm^-3. PTF (Rawls) substantially overestimated both the field capacity(0.413 cm³ cm^-3) and the permanent wilting point (0.176 cm³ cm^-3), resulting in an overestimation of +33% for the plant-available water capacity. PTF (Saxton) performed well on the field capacity(0.302 cm³ cm^-3) but underestimated the permanent wilting point (0.115 cm³ cm^-3), resulting in a 5% overestimation of plant-available water capacity. The sensor approach underestimated the field capacity (0.204 cm³ cm^-3) and remained close to the permanent wilting point (0.127 cm³ cm^-3), resulting in a 56% underestimation of plant-available water capacity. At a depth of 50 cm, the Lab-measured volumetric water content at the field capacity was 0.238 cm³ cm^-3, and at the permanent wilting point, it was 0.174 cm³ cm^-3, resulting in a plant-available water capacity of 0.064 cm³ cm^-3. PTF (Rawls) overestimated the field capacity (0.357 cm³ cm^-3), while understating the permanent wilting point (0.156 cm³ cm^-3), leading to a plant-available water capacity of 0.201 cm³ cm^-3. Similar tendencies at 50 cm can be observed for PTF (Saxton), resulting in a plant-available water capacity of 0.185 cm³ cm^-3. The Sensor measurements at a depth of 50 cm overestimated the field capacity (0.282 cm³ cm^-3), resulting in a plant-available water capacity of 0.038 cm³ cm^-3. At a depth of 130 cm, the field capacity was 0.334 cm³ cm^-3, and the permanent wilting point was 0.304 cm³ cm^-3, yielding a low plant-available water capacity of 0.03 cm³ cm^-3. PTF (Rawls) and PTF (Saxton) significantly overestimated the plant-available water capacity at 130 cm (>0.170 cm³ cm^-3), whereas the Sensor approach underestimated the field capacity (0.272 cm³ cm^-3) and the permanent wilting point (0.250 cm³ cm^-3), resulting in a low plant-available water capacity of 2.7 cm³ cm^-3. The PTF (Rawls) overestimated the plant-available water capacity throughout the whole soil profile. With the absence of organic carbon, at 50 cm downwards, the PTF (Saxton) aligned with the PTF (Rawls). The Sensor approach, on the other hand, underestimated the plant-available water capacity throughout the soil profile.

Figure 5 illustrates the impact of the three approaches used to generate the field capacity and permanent wilting point (Lab-measured, PTF Saxton, and PTF Rawls) on the simulated soil water content. The simulations are compared with observed values at a depth of 15 cm, and model accuracy is evaluated using RMSE for the years 2022, 2023, and 2025. In 2022, all three models closely tracked the observed gradual decline in soil moisture from approximately 0.30 to 0.15 cm³ cm^-3 between May and July, with RMSE values lower than 0.037 cm³ cm^-3 across all approaches. A significant best fit was observed for the approach PTF (Rawls) (0.011 cm³ cm^-3). The 2023 season was characterized by higher moisture variability, with spring and summer peaks, which increased RMSE across all approaches. The lowest RMSE of 0.0363 cm³ cm^-3 was observed for the PTF (Rawls), followed by the Lab-Measured (RMSE of 0.0467 cm³ cm^-3) and PTF (Saxton) (RMSE of 0.0604 cm³ cm^-3). In 2025, the PTF (Rawls) showed the lowest RMSE of 0.0266 cm³ cm^-3, whereas the RMSE for the PTF (Saxton) and the Lab-Measured were>0.06 cm³ cm^-3. Overall, PTF (Rawls) provided the most consistent simulation accuracy across all three seasons, and greater moisture variability reduced model performance. Figure 6 compares observed soil water content with simulated values from the three approaches: Lab-Measured, PTF (Saxton), and PTF (Rawls) at a depth of 50 cm across the 2022 and 2023 growing seasons. The observed values consistently remain higher than all simulated outputs throughout both years, suggesting a general tendency of the models to underestimate soil moisture. In both years, soil water content shows a clear declining trend from spring through summer, reflecting seasonal drying of the soil profile. The PTF (Rawls) approach achieved the lowest RMSE in both years (0.053 cm³ cm^-3 in 2022 and 0.045 cm³ cm^-3 in 2023). The Lab-Measured approach showed intermediate performance with RMSEs of 0.068 cm³ cm^-3 in 2022 and 0.043 cm³ cm^-3 in 2023, while PTF (Saxton) had the highest errors at 0.078 and 0.07 cm³ cm^-3, respectively. Overall model performance improved slightly from 2022 to 2023 across all three approaches, but performance at 50 cm was weaker than at 15 cm, as indicated by lower RMSE.

3.3 Evaluation of model performance of the soil profile optimization

In Table 3, model performance was assessed across four simulation approaches: PTF (Rawls), PTF (Saxton), Bayesian, and Grid search. For the inverse optimization algorithms, the approaches are further distinguished by their specific calibration objective (soil water content, above-ground biomass, or grain weight) and subsequently evaluated across all individual validation metrics and overall model performance. The default PTF (Rawls) approach showed the lowest simulation errors across all three variables and the best model performance (nRMSE of 5%), while the PTF (Saxton) showed the highest errors across all metrics (nRMSE of 20%). Compared to the PTF (Saxton), all inverse optimization approaches improved model performance. Among these optimization approaches, the model performance (nRMSE) ranged from 16% (Bayesian targeting above-ground biomass) to 9% (Grid search targeting above-ground biomass). When evaluating variable-specific validation, the highest simulation errors for soil water content were observed when the models were optimized against other crop metrics (nRMSE of 23.3% for Bayesian targeting above-ground biomass, and 11.5% for Grid search targeting grain weight). Conversely, the best optimization approach for simulating soil water content at 15 cm was the Grid search targeting soil water (nRMSE of 3.6%). For above-ground biomass validation, the Bayesian approach targeting above-ground biomass showed the lowest error among the optimizers (nRMSE of 10.4%). Finally, grain weight was best simulated by the Bayesian approach targeting grain weight (nRMSE of 7.4%), followed closely by the Grid search targeting above-ground biomass (nRMSE of 7.7%). Ultimately, the best overall optimization model performance across all variables was achieved by the Grid search targeting above-ground biomass (nRMSE of 9%).

Table 3. Evaluation of model performance across four simulation approaches (PTF (Rawls), PTF (Saxton), Bayesian, and grid search).

For the inverse modeling approaches, separate calibration runs were performed targeting individual variables (Target variable). The resulting models were evaluated across all observed metrics (Validation Metrics) over 2022, 2023, and 2025.

The two optimization approaches —Grid search (targeting above-ground biomass) and Bayesian (targeting grain weight)—were selected because they demonstrated the lowest simulation error when targeting these crop growth parameters. In Table 4, these approaches are used to demonstrate the impact of the optimization procedure on the resulting van Genuchten and functional soil hydraulic parameters. The results of the optimization approaches are compared with both the direct lab-measured data and the PTF (Saxton) estimations, which serve as the baseline scenario and the initial starting point for the optimization process. Throughout the optimization process, the structural van Genuchten shape parameters (α and n) as well as the saturated and residual water contents (θ_r and θ_s) remained rigid. Both the Grid search (above-ground biomass) and Bayesian (grain) optimization algorithms resulted in values identical to the PTF (Saxton) baseline across all horizons. This indicates that the curve-fitting procedure could not mathematically differentiate a better parameter combination, even though the resulting R² remained consistently high. However, when interpreting this lack of deviation, it must be considered that only three measured values were available for the fitting process, which inherently limited the algorithm’s resolution to further optimize these specific structural parameters.

Table 4. Comparison of the simulation approaches PTF (Saxton), grid search targeting above-ground biomass, and Bayesian targeting grain weight with the lab-measured values for the derived van Genuchten parameters, their evaluation (R²), and the soil hydraulic parameters: field capacity and permanent wilting point and plant-available water.

While the fitted van Genuchten parameters remained unchanged, the optimization approaches induced shifts in the hydraulic parameters. Both optimization algorithms increased the field capacity across all five horizons compared to the initial PTF (Saxton) and the Lab-Measured values. A divergence emerged between the two optimization approaches regarding the permanent wilting point. The Grid search (above-ground biomass) maintained a permanent wilting point identical to the PTF Saxton baseline throughout the entire soil profile. In contrast, the Bayesian (grain) introduced horizon-specific variations, reducing the permanent wilting point in the Ap1, Big1, and Big2 layers, resulting in plant-available water being higher in the optimized profiles. The Bayesian (grain) increased the plant-available capacities, driven by the combined effect of elevated field capacities and selectively lowered permanent wilting points (0.250 cm³ cm^-3 in the Ap1 and 0.263 cm³ cm^-3 in the Big2 horizon). The Lab-measured data exhibited significant discrepancies from both the PTF (Saxton) and the optimized estimations, particularly in the deeper subsoil layers.

3.4 Impact of simulation approaches on output parameters of the soil water balance

Table 5 compares the output parameters from the soil water balance, including cumulative transpiration, soil evaporation, evapotranspiration, and extractable soil water (in mm), across four approaches: Lab-Measured, PTF (Rawls), Bayesian (Grain), and Grid search (above-ground biomass). While cumulative precipitation remained relatively consistent (ranging from 407 to 510 mm), the Lab-Measured approach consistently yielded the lowest cumulative transpiration and evapotranspiration values. The cumulative transpiration for PTF (Rawls) was consistently lower than both optimization approaches, Bayesian (grain) and Grid search (above-ground biomass), where the 2025 season showed a higher cumulative transpiration of 68% for the Bayesian (grain) and 52% for the Grid search (above-ground biomass) compared to PTF (Rawls). The 2023 season showed the highest cumulative transpiration, with a mean of 264 mm across all simulation approaches. For cumulative soil evaporation, no significant differences were observed between the simulation approaches. Similar to the cumulative transpiration, the seasonal trend with the highest evaporation was simulated for 2023. Cumulative evapotranspiration showed a trend similar to that observed for transpiration, with the highest evapotranspiration found for Bayesian (Grain)>Grid search (Above-ground biomass)>PTF (Rawls)>Lab-Measured. The mean percentage increase in the cumulative evapotranspiration compared to PTF (Rawls) was 12% for Bayesian (Grain) and 7% for Grid search (Above-ground biomass). The extractable soil water at a 50 cm depth was evaluated at the onset of water stress, defined as the first day of the growing season when pF exceeded 2.7. Across simulation approaches, the Lab-Measured values showed the lowest extractable soil water content, at 19 mm at 221 days after planting (DAP) in 2022 and 42 mm at 224 DAP in 2023. The Bayesian (Grain) approach simulated the highest levels of extractable soil water, while the PTF (Rawls) and Grid search(Above-ground biomass) approaches showed intermediate soil water contents.

Table 5. Comparison of simulated soil water balance output parameters (cumulative precipitation, transpiration, soil evaporation, evapotranspiration and extractable soil water (mm) at pF>2.7 at 50 cm depth) across four simulation approaches: lab-measured, PTF (rawls), Bayesian (grain), and grid search (above-ground biomass) at specific days after planting (DAP) for the years 2022, 2023, and 2025.

Figure 7 illustrates the simulated daily water stress factor for transpiration across the 2022, 2023, and 2025 growing seasons, comparing four distinct simulation approaches, with 0 denoting no water stress and 1.0 indicating maximum water stress. Across all three years, the Lab-Measured approach is the most sensitive indicator, consistently predicting the onset of water stress between mid-May and early June. In contrast, the three alternative simulation methods—PTF (Rawls), Bayesian (Grain), and Grid search (Above-ground biomass)—exhibit the stress response later in the season. The in-season matric potential measurements detected water stress in 2022 at 221 DAP and in 2023 at 224 DAP. Compared to that, the Lab-Measured approach showed the first peak in water stress at 186 DAP(35 days earlier), PTF (Rawls) at 220 DAP(1 day earlier), Grid search at 232 DAP (11 days later), and Bayesian (Grain) at 240 DAP (19 days later). In 2023, the Lab-Measured approach simulated water stress at 194 DAP (30 days earlier), PTF (Rawls) at 213 DAP (9 days earlier), Bayesian (Grain) at 218 DAP (6 days earlier), and Grid search(Above-ground biomass) at 234 DAP (10 days later).

4 Discussion

The application of the DSSAT soil water model relies on the tipping bucket method, which addresses well-documented structural limitations inherent in empirical modeling. While critics often point to its simplified physics, the primary advantage of the DSSAT soil module lies in its operational 'easy-to-use' architecture, which prioritizes computational robustness and low data requirements—making it suitable for multiple modeling scenarios. Unlike the mechanistic models that solve the Richard Equation, the tipping bucket approach avoids numerical instability and convergence issues, providing a stable platform for synchronizing soil water status with crop growth modules across diverse environments (Farthing and Ogden, 2017). Despite these operational strengths, the challenges associated with the tipping bucket approach include the use of PTF, the runoff curve number method, and water uptake mechanisms (Camargo and Kemanian, 2016; Gijsman et al., 2002; Sadler et al., 2000). As stated by Wallach et al. (2024), reducing simulation errors in such models typically follows one of two trajectories: either a fundamental redesign of model structure or a rigorous refinement of parametrization procedures. Jarvis et al. (2022) questioned the effectiveness of a tipping-bucket approach in crop models because it cannot simulate upward water movement, does not account for evaporation at the soil surface and for capillary rise into the root zone, and thus leads to inaccurate estimates of crop transpiration during dry summer periods. However, adjusting and changing these kinds of model algorithms requires a depth of data that is rarely available. The dataset for this field study, with its wheat growth data and daily measurements of volumetric soil water content and matric potential at multiple depths, is therefore unique and rare. Soil matric potential is considered the superior indicator for detecting plant water stress because it directly measures the energy state of soil water, representing the force roots must overcome to extract water, and its standard values are readily available for individual crop species and soil textures (Campbell, 1985; Taylor, 1965). Conversely, soil water content remains a critical variable for modeling and monitoring, as it enables continuous measurement of water dynamics throughout the root zone, but lacks a direct link to the plant’s current status. Despite its utility, controlling the upper and permanent wilting points based on volumetric soil water content measurements is complex, as seen in the current data set, where even under drought conditions, volumetric water content did not decrease. The water stress bands in Figures 2, 3 were adapted from Feddes et al. (1978), which relates root water uptake and transpiration demand to matric potential, specifically defining growth reduction between pF 2.7 and 3.0. However, these bands remain a proxy, as they are inherently dependent on cultivar-specific transpiration rates. To improve future model structures, data on root growth dynamics and specific crop water potential should be incorporated; such data collection is resource-intensive. Consequently, this study focused on parameterizing soil profile characteristics and demonstrated a potential ready-to-use strategy through a systematic optimization procedure.

A fundamental challenge in environmental modeling is that soil processes do not respect national boundaries, yet the data describing them remain fragmented by local standards. In this study, determining soil texture classes proved challenging due to discrepancies between the national classification and international standards. By comparing the German standard ISO 11277:2020 and the USDA standard, the shift is evident in the silt texture class: both specify a minimum size of 2 μm, but the maximum size varies (e.g., 63 μm in ISO 11277:2020 and 50 μm in USDA). While these definitions are well-documented in the literature, they are frequently dismissed as minor technicalities (Nemes et al., 1999). This shift is a problem that often escapes modelers’ notice, resulting in systemic bias in hydraulic parametrization. Especially the soil used in this study, characterized by a high silt content, leads to a significant shift by adapting from ISO 11277 to USDA, specifically, a shift towards higher sand content of ≥6% in the topsoil (0–18 cm) and 4% in the subsoil (18–69 cm). The adaptation was performed using the Moyes approach, which is a proxy for recalculating the texture fractions. However, as noted in the development of the HYPRES database, mathematical interpolation is often a necessary second-best to direct measurements when reconciling disparate datasets (Nemes et al., 1999). The significance of this textural shift becomes critical when these data serve as primary inputs for PTFs. Since PTFs are statistical models that translate basic soil properties into hydraulic parameters, any systematic bias in the input data propagates directly into the predicted soil hydraulic parameters. Nemes and Rawls (2006) addressed this issue, demonstrating that harmonizing size boundaries can lead to substantial errors in predicted retention curves. They suggest that using interpolated data—despite the inherent magnitude of interpolation error—poses fewer risks when applying a PTF than using measured data that follow inconsistent silt-to-sand boundaries.

When comparing the retention curves derived from the default PTFs of Saxton and Rawls with the curves fitted to the lab- and sensor-measured data points, significant differences were detected across the entire soil profile, as evidenced by the PTF-based approaches overestimating the Lab-Measured values. In their foundational work on soil water availability, Ratliff et al. (1983) argued that laboratory-measured field capacity is generally inadequate for accurate soil-water modeling because it rarely matches in situ field conditions. In theory, the field capacity at -33 kPa is used for moderately coarse and finer textured soils, whereas -10 kPa is used for coarse-textured soils (Jamison and Kroth, 1958). To date, little literature has examined how crop growth models address the impact of hydraulic parameters. Yao et al. (2020) used the DSSAT-CSM-CERES-Wheat and defined the field capacity at -10 kPa for a silt loam. However, in this study, following the definitions of Rawls and Saxton, a matric potential of -33 kPa (pF 2.5) was defined as the field capacity, which is also consistent with the German KA5 logic, as it characterizes an unaffected soil by groundwater. Compared with the retention curves from the Lab-Measured data, both PTF (Saxton) and PTF (Rawls) deviated significantly, with a strong tendency toward overestimation. While Patil and Singh (2016) reported RMSEs for PTFs ranging from 0.007 to 0.07 m³ m^-3, the error values in the current study exceeded these benchmarks for estimated hydraulic parameters. Especially for soil texture classes with high silt content, precise determination of soil water content at the field capacityof -33 kPa is reported to be prone to overestimation, which may also be attributed to the original datasets lacking silt-dominated texture classes. Plant-available water capacity is converted to extractable soil water (plant-available water) by scaling it to the soil profile depth and the crop’s rooting zone. Calculating the extractable water via the mean values defined for silt loam (0–50 cm), silty clay loam (50–69 cm), and silty clay (69–100 cm) after Allen et al. (1998) would result in a total of 135 mm m^-1, whereas the lab-measured values resulted in a total of 92 mm m^-1, indicating that the current soil generally showed low values compared to standards reported in literature.

In the current study, significant challenges in simulating soil water as a state variable were addressed. In these simulations, several important tendencies for future model runs were observed. The integration of C_org with the PTF (Rawls) resulted in the highest predictive accuracy for simulating soil water content at a depth of 15 cm. The selection of the Rawls PTF is automatically triggered within the model when C_org is used as an input. Rawls et al. (2003) stated that incorporating C_org into their regression model improved performance by 25% compared to using textural classes alone. In this study, the PTF (Rawls) demonstrated high simulation accuracy for soil water content at 15 cm, achieving an nRMSE of 2.9%—a 12% improvement over the PTF (Saxton). These results suggest that incorporating C_org measurements is recommended to leverage the PTF’s (Rawls) superior accuracy for precise topsoil water simulations. If C_org is not available, the inverse modeling approach using both optimization approaches, Bayesian and Grid search, showed promising results, significantly improving overall model performance by targeting each of the target variables. Merely targeting the above-ground biomass with the Bayesian algorithm resulted in higher nRMSE for the soil water content. Generally, selecting suitable target variables and setting the optimization range for the selected hydraulic parameters are complex tasks. While the Bayesian Optimization was constrained to nominal bounds of 15% and the Grid Search to 20%, this distinction accounts for the structural differences between continuous and discrete algorithms. Because Bayesian optimization explores parameters independently within a continuous volume, opposing shifts of 15% in the field capacity and permanent wilting points can functionally alter the plant-available water capacity by up to 30%. Grid search, which is restricted by discrete intervals and constrained combinations, requires slightly wider nominal boundaries to encompass a comparable functional parameter space (Bergstra and Bengio, 2012). This structural efficiency translates into a significant difference in computational cost. For the current study, the Grid search required approximately three hours for the optimization procedure, whereas the Bayesian optimization converged in under ten minutes. When extrapolating this computational demand to site-specific field management—where capturing spatial soil variability requires independent, repeated calibrations across numerous management zones or grid cells—the exponential computational burden of grid search quickly becomes prohibitive. Thus, the accelerated convergence of Bayesian optimization is not merely a methodological preference, but a practical necessity for scalable, high-resolution agricultural modeling.

Overall, the soil-water simulations for the 50 cm soil layer showed higher errors than those for the topsoil simulations. This is consistent with the observations of Braga and Jones (2004), who also reported underestimation of soil water content at greater depths. In this case, it can be assumed that the tipping bucket failed to capture upward flux, a phenomenon that is particularly pronounced in silty soils because their pore sizes are small enough to create high suction while still allowing water to move at a reasonable speed. Optimizing the soil water content targeting 50 cm did not improve the model, indicating that the model dynamics cannot be captured by the tipping bucket (Jarvis et al., 2022). The impact of the different approaches on the output parameters that drive the soil water balance showed that differences in evapotranspiration were driven by plant transpiration rather than by soil evaporation. Furthermore, simulating the plant’s water-stress response is not driven by the soil; it is driven by the wheat-specific stress response. This was evident in the comparison with the lab-measured treatment, which showed a strong water-stress response in all three years. The simulated extractable water at the onset of water stress was significantly lower for the lab-measured approach in both years than for the other treatments. The model demonstrated high structural ability, maintaining consistent outputs even under rigorous optimization. The current absolute values deviate from the observed stress response; the trend remains internally consistent. However, these results showed general limitations in the model’s simulation of crop responses to water stress, indicating that upward movement is crucial for controlling evaporation and for the capillary rise of groundwater into the root zone, which supports crop transpiration (Longo et al., 2021; MacBean et al., 2020). While the results in this study highlight the inherent challenges in modeling soil hydraulic parameters, they underscore that moving beyond legacy empirical frameworks toward more mechanistic soil models is a priority for capturing the soil-crop interface with the required precision.

5 Conclusion

This research provides a framework for reducing predictive uncertainty in crop models by evaluating the standard PTFs in DSSAT using a detailed experimental dataset. By differentiating between the use of PTF Saxton and PTF Rawls, the advantage of measured C_org was evident, as indicated by the best goodness of fit, with an nRMSE of merely 5% over three years of simulations. However, in the absence of C_org data, model simulations rely on the Saxton approach, resulting in greater inaccuracy and an nRMSE of 20%. To address the uncertainty, a ready-to-use optimization approach was implemented via the SPO software plug-in using inverse modeling to optimize selected hydraulic parameters by targeting either crop growth or water content. The SPO incorporates Bayesian and Grid Search algorithms, significantly improved for both targeting growth and soil data; specifically, Bayesian optimization achieved an overall model performance of 10% nRMSE for grain weight, while the Grid Search approach performed best for above-ground biomass with an nRMSE of 9%. Despite these advancements, the study identified critical boundaries in simulating subsoil water content. Specifically in silt-dominant soils, uncaptured upward flow dynamics notably affected soil water predictions, indicating that future model updates warrant the inclusion of a fully mechanistic soil water module. Finally, to establish the broader applicability, future studies must evaluate it across diverse soil types and textural classes, as differing pedological and structural conditions will dictate distinct parameterization strategies.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. During the preparation of this work, the author used Google Gemini to troubleshoot data visualization in R, code in Python, and Grammarly to improve the manuscript’s readability and language. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the published article.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fagro.2026.1912611/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Tiago B. Ramos, National Institute for Agricultural and Veterinary Research (INIAV), Portugal

Reviewed by: Raúl Sánchez, Polytechnic University of Madrid, Spain

Menglu Huang, The Hong Kong Polytechnic University, Hong Kong SAR, China