# Inverse Modelling to Estimate Soil Hydraulic Properties at the Field Scale

Authors: Mohamed Khaled Salahou; Xiaoyuan Chen; Yupeng Zhang; Xiyun Jiao; Haishen Lü
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Field-scale estimation of soil hydraulic parameters is important for simulating water movement, particularly in the vadose zone. The accuracy of soil hydraulic parameters significantly affects the output of models using these parameters as input variables. A new approach for an inverse method, often referred to as inverse modelling, was presented to calibrate soil hydraulic parameters. The new method is based on the use of uniform design theory to generate uncertainty estimates of the soil hydraulic parameters. The uncalibrated and calibrated soil hydraulic parameters were evaluated by comparing the predicted and measured soil water contents using the HYDRUS-1D model. The results show that the calibrated soil hydraulic parameters that were inversely estimated by uniform design theory were accurate and reasonable. In addition, the results provided recommendations for improving the estimation of soil hydraulic parameters when there are no available observation data for calibrating and optimizing soil hydraulic parameters.


---

## Body

## 1. Introduction

The accurate estimation of soil hydraulic parameters is crucial in simulating water movement [1–4], as well as in providing input for related agro-hydrological models (i.e., the Green-Ampt model [5–7]). These simulation models require accurate information about the unsaturated soil hydraulic conductivity and the soil water retention curve. The precision of estimations of soil hydraulic parameters has increased considerably during the last decade by improving the methods used to predict the required soil hydraulic parameters from the soil properties [8–10], as direct measurements of soil hydraulic characteristics are costly and time-consuming. Relationships between soil hydraulic parameters and other (textural) properties are commonly called pedotransfer functions (PTFs) [11].PTFs have been developed using various mathematical and statistical approaches, such as neural network analyses or regression, and these functions can be used to directly predict soil hydraulic characteristics or to predict the parameters used in analytical models [12–16] to determine soil hydraulic characteristics. The Brooks and Corey (BC) and the (VG) models are commonly used analytical models [17–22]. The BC parameters are presented in Rawls et al. [23] who used multiple linear regression for a large database of 2540 soil horizons. The VG parameters are predicted from ROSETTA, which is a software program that uses hierarchical neural network analysis to predict VG soil hydraulic parameters; an increase in accuracy is obtained when additional input variables are used. ROSETTA implements five PTFs for the estimation of soil hydraulic properties using limited (textural class) to more extensive (texture, bulk density, and one or two water retention points) information [24]. The advantage of using neural networks to establish PTFs is that no a priori concept is needed. In addition, ROSETTA is very user friendly [25].In contrast, inverse modelling aims to calibrate soil hydraulic parameters by minimizing the difference between observed and simulated data, such as the metric potential [26], water soil content [22, 27–30], solute concentration, and the combination of in situ observation variables [31, 32]. Previous studies that have used different forms of inverse modelling to calibrate soil hydraulic parameters include the Gauss-Marquardt-Levenberg method [33] (i.e., [28, 30]), the Shuffled Complex Evolution Metropolis algorithm (SCEM) [34] (i.e., [35, 36]), and NSGA-II [37] (i.e., [29]). Most of these models are based on Monte Carlo sampling. Uniform design (UD) theory, which was proposed by [38–41], is also based on the quasi-Monte Carlo method or number theoretic method. The UD method was developed to select a small but representative set of points from a given set of points such that the sampled points are uniformly scattered. The UD method has attracted attention in many studies in various fields [42–46] because it has the following distinct characteristics: space filling, robustness, and multiple levels [40, 41]. These characteristics have led to the UD method becoming one of the major current experimental design approaches in China. Additionally, the UD method offers many experimental tables for users to conveniently utilize. For these reasons, we would like to introduce the UD method as a new parameter sampling method. Thus, the objectives of this study were as follows: (1) to determine soil hydraulic parameters using different soil physical information (uncalibrated parameters); (2) to calibrate and optimize the soil hydraulic parameters based on a comparison of predicted and measured soil water contents using the HYDRUS-1D model, where the uncertainty estimates were generated by UD theory; and (3) to determine which levels of soil physical information are sufficient for estimating soil hydraulic parameters. (The calibrated and uncalibrated soil hydraulic parameters were used to compare the predicted and measured soil water contents using the HYDRUS-1D model).


## 2. Methods

### 2.1. Study Area

The field experiment was conducted at the CAS Ecological Agricultural Experiment Station of Hebei, China (116°40′ E, 38°06′ N). This area is in the monsoon climate zone. The average temperature and altitude are approximately 12.3°C and 20 m, respectively. The mean annual precipitation is approximately 567.4 mm, and annual evaporation ranges between 1500 and 1800 mm. Rainfall generally occurs during the summer, with 73%, 11%, 13%, and 3% of the annual precipitation occurring during the summer, spring, autumn, and winter, respectively. The site has a silt loam soil and a groundwater depth of approximately 4 m. The soil physical characteristics, including bulk density and particle size distribution, were determined in the laboratory for representative depths, as shown in Table1.Table 1
Physical soil properties at the experimental site based on the USDA classification system.
Soil depth (cm)Particle size distribution (%)Soil textureBulk density (g.cm−3)SandSiltClay0–6011.8466.2121.95Silt loam1.4060–27517.4465.7016.86Silt loam1.42275–3456.8660.9532.19Silty clay loam1.44345–44565.0927.657.26Sandy loam1.44445–65085.5410.434.03Sandy soil1.20


### 2.2. Model Simulation Using HYDRUS-1D

Five different scenarios were assumed to obtain the soil hydraulic parameters. The HYDRUS-1D model has a database for soil hydraulic parameter estimation (uncalibrated parameters); these data were originally obtained from Rawls et al. [23] for the BC model and from ROSETTA software [47] for the VG model. Thus, the soil hydraulic parameters for the first four scenarios (S1, S2, S3, and S4) were selected from the HYDRUS-1D database of soil hydraulic properties based on the hydraulic model (BC or VG) and different soil physical information levels. To compare the uncalibrated and calibrated soil hydraulic parameters, the measured profile soil water content was used to calibrate and optimize the soil hydraulic parameters by inverse modelling (scenario S5). Finally, the obtained parameters from each previous scenario (calibrated and uncalibrated parameters) were fed into the HYDRUS-1D model to determine the accuracy of each scenario parameter by comparing the observed and simulated soil water contents. The details of the VG, BC, and other HYDRUS-1D model parameters are as follows.


#### 2.2.1. Water Flow Equation

The transport of soil water was numerically simulated by solving the Richards equation (48) using the model developed by the United States Salinity Laboratory (HYDRUS-1D) [49, 50]. The HYDRUS-1D model solves the following one-dimensional Richards equation for nonsteady unsaturated water flow:(1)∂θ∂t=∂∂zKh∂h∂z+1−S,where h is the soil water pressure head [L], θ is the volumetric soil water content [L3.L−3], t is time (T), z is the soil depth [L], K(h) is the unsaturated hydraulic conductivity function [L.T−1], and S is the sink term accounting for root water uptake [L3.L−3.T−1].


#### 2.2.2. Soil Hydraulic Model and Properties

To solve the Richards equation, the water retention functionθ(h) and the hydraulic conductivity function K(h) were obtained from the BC and VG hydraulic models.(1) Brooks and Corey Model. The unsaturated soil hydraulic properties θ(h) and K(h) in (1), according to Brooks and Corey [12], are given by(2)θh=θr+θs−θrαh−nh<−1α,θsh≥−1α,Kh=KsSe2/n+l+2h<−1α,Ksh≥−1α,Se=θh−θrθs−θr,where Se is the effective saturation; θr is the residual moisture content, expressed as [L3.L−3]; θh,θs,Ks, and h are as previously defined; α is the inverse of the air entry point; n is the pore size distribution index; and the HYDRUS-1D database of soil properties provides a value of l=1 for all BC soil data sets. The parameters α, n, and l are considered to be empirical coefficients that affect the shape of the hydraulic functions. As mentioned above, the BC parameters were originally reported by Rawls et al. [23].(2) Van Genuchten-Mualem Model. The unsaturated soil hydraulic properties are described by the VG model [15, 51]:(3)θh=θr+θs−θr1+αhnmh<0,θsh≥0,Kh=KsSel1−1−Se1/mm2h<0,Ksh≥0,where m=1−1/n,n>1 [52], and l was estimated by [51] as an average value of approximately 0.5 for many soil types; other variables are as previously defined. The VG parameters were obtained using ROSETTA software (VG-ROSETTA) with input data as follows: soil textural class; particle size distribution; or particle size distribution and bulk density.


#### 2.2.3. Evapotranspiration

Meteorological data, such as air temperature, solar radiation, wind speed, precipitation, and relative humidity, were collected from the Changzhuang measuring station, located near the experimental site (116°51′ E, 38°21′ N), for the calibration period (4th June 2009 to 12th November 2009) and the validation period (7th May 2010 to 1st October 2010). The daily reference evapotranspiration ET0 (L.T−1) was calculated using the FAO-56 Penman–Monteith equation (53). The potential evapotranspiration rate, ETPt, was calculated by multiplying the crop coefficient (Kc) by the ET0 rate and was divided into the potential evaporation rate, EPt, and potential transpiration rate, TPt, by Beer’s Law (Ritchie [54]:(4)EPt=ETPte−σLAI,TPt=ETPt−EPt,where LAI is the leaf area index with units of [L2.L−2], and σ is a plant-dependent parameter (dimensionless), for which Ritchie [54] suggested a value of 0.4 for agricultural crops with an LAI less than 2.7 cm2 cm−2. The Kc and LAI values for a cotton crop were obtained from [55].


#### 2.2.4. Root Water Uptake

The Feddes model was used to estimate the water uptake function of the root system [56].(5)Sh,z=γhβzTPt,where S(h, z) is the sink; TPt is the potential transpiration rate [L.T−1]; βz is the root density distribution function (L.L−2) (the integration of the root distribution along the root zone is equal to one; there are many ways to express the βz function, and HYDRUS-1D does not limit its users in any way); and γh is the soil water stress response function (dimensionless) of Feddes et al. [56].(6)γh=0;h≥h1,h1−hh1−h2h1>h≥h2,1;h2>h≥h3,h−h4h3−h4h3>h≥h4,0;h4≥h,.The values of the threshold parameters (h1, h2, h3, and h4) in the cotton root water uptake model were set to h1 = −1 cm, h2 = −25 cm, h3max = −200 cm, h3min = −800 cm, and h4 = −8000 cm. These values were optimized by Li [57] for cotton water uptake simulation.


### 2.3. Inverse Model Based on Uniform Design Theory

The HYDRUS-1D model has a database for soil hydraulic property estimation, and the estimated hydraulic parameters of BC and VG-ROSETTA should be accompanied by uncertainty estimates that permit an assessment of the reliability of parameter predictions. UD theory [39, 58] was used to generate the uncertainty estimates of the soil hydraulic parameters. Fang [39] made use of so-called good lattice point (glp) sets to construct many tables for UD. UD tables are denoted as follows:(7)Unqs,where U denotes the UD, n is the number of runs, s is the number of factors, and q is the number of levels.A UD table design by Fang [39] with 2 < s < 30, 3 < q < 11, and different numbers of runs (n) was used. For a large number of factors (s > 10), the number of runs is n = 110, 121, 132, 143, 154, or 165. The UD table can also be obtained from the Hong Kong Baptist University Mathematical Department website (https://www.math.hkbu.edu.hk/UniformDesign/).An uncertainty level of ±25% of the initial value of each soil hydraulic parameter was analysed; hence, each parameter was set at 11 levels (q = 11) (i.e., αinitial,αinitial±0.05,αinitial±0.10,αinitial±0.15,αinitial±0.20,andαinitial±0.25), and there was a total of five soil layers (Table 1), each of which consisted of five parameters (α, n, θs, θr, and Ks). Therefore, a total of s=5×5=25 parameters needed to be estimated. A UD table of 110 runs with 25 factors and 11 levels, designed by Fang [39], was selected U1101125. The UD test tables were obtained through the conversion of the level and actual value of each parameter. Then, the 110 sets of parameters were fed into the HYDRUS-1D model. The parameters that resulted in the highest goodness-of-fit between the observed and simulated soil water content were selected.The soil water contents were measured by time-domain reflectometry (TDR) probes on a daily basis during the calibration period (4th June 2009 to 12th November 2009) and the validation period (7th May 2010 to 1st October 2010) in soil layers at depths of 50, 100, 150, 200, 250, 300, 350, 400, and 450 cm in the study area.Table2 summarizes the input and output parameters for each scenario. For scenario S1, where the listed source is [23], parameters were created for soil texture. For scenarios where the listed source is VG-ROSETTA, parameters were derived from ROSETTA [47] using the soil textural class option for scenario S2, the particle size distribution option (sand, silt, and clay fractions) for scenario S3, and the particle size distribution and bulk density option for scenario S4. For scenario S5, the initial values of the parameters α,n,θs,θr,andKs were obtained from ROSETTA using the particle size distribution and bulk density as input data; subsequently, an uncertainty level of ±25% of the initial value of each parameter was analysed; the parameters from the UD table U1101125 that resulted in the highest goodness-of-fit between the observed and simulated soil water content were selected as the output parameters.Table 2
Input and output parameters of the model options.
ScenarioHydraulic modelSourceInput∗Output parametersParameter typeS1BC[23]Textural class (USDA)α, n, θs, θrKsUncalibratedS2VGVG-ROSETTA∗∗Textural class (USDA)α, n, θs, θrKsS3VGVG-ROSETTASand, silt, and clay fractionsα, n, θs, θrKsS4VGVG-ROSETTASand, silt, and clay fractions + soil bulk densityα, n, θs, θrKsS5VG——α, n, θs, θrKsCalibrated∗ input data for soil physical characteristics were obtained from Table 1. ∗∗ For scenarios where the listed source is VG-ROSETTA, parameters were derived from ROSETTA software [47].


### 2.4. Initial and Boundary Conditions

The initial condition was set as the initial soil water content in the study domain at timet = 0. The upper boundary condition represents the atmospheric boundary condition, where water flux depends on precipitation, evaporation, and the moisture conditions in the soil. The lower boundary condition was the measured daily pressure head derived from the groundwater table depth. The position of the groundwater table was measured every day during the calibration and validation period in an observation well located inside the experimental field (Figure 1). The ETPt values, groundwater depths, and precipitation rates during the calibration and validation periods are shown in Figure 1.Figure 1
ETPt, groundwater depth, and precipitation rate in (a) 2009 and (b) 2010.
(a)(b)


### 2.5. Statistical Analysis of the Results

Two statistical indictors were used to analyse the goodness-of-fit between the observations and simulations: (1) the root mean square error (RMSE) and (2) the Nash-Sutcliffe modelling efficiency (NSE) [59], which are calculated as follows:(8)RMSE=1N∑i=1NPi−Oi2,NSE=1−∑i=1NPi−Oi2∑i=1NOi−O¯2,where Pi represents the predicted value, Oi represents the observed value, O¯ is the mean value of all observed events at time i, and N is the number of data pairs. The RMSE provides an indication of how good the agreement between the predicted and observed values is. NSE ranged between -∞ and 1, where a high value of NSE (NSE > 0) represents a good model performance.


## 2.1. Study Area

The field experiment was conducted at the CAS Ecological Agricultural Experiment Station of Hebei, China (116°40′ E, 38°06′ N). This area is in the monsoon climate zone. The average temperature and altitude are approximately 12.3°C and 20 m, respectively. The mean annual precipitation is approximately 567.4 mm, and annual evaporation ranges between 1500 and 1800 mm. Rainfall generally occurs during the summer, with 73%, 11%, 13%, and 3% of the annual precipitation occurring during the summer, spring, autumn, and winter, respectively. The site has a silt loam soil and a groundwater depth of approximately 4 m. The soil physical characteristics, including bulk density and particle size distribution, were determined in the laboratory for representative depths, as shown in Table1.Table 1
Physical soil properties at the experimental site based on the USDA classification system.
Soil depth (cm)Particle size distribution (%)Soil textureBulk density (g.cm−3)SandSiltClay0–6011.8466.2121.95Silt loam1.4060–27517.4465.7016.86Silt loam1.42275–3456.8660.9532.19Silty clay loam1.44345–44565.0927.657.26Sandy loam1.44445–65085.5410.434.03Sandy soil1.20


## 2.2. Model Simulation Using HYDRUS-1D

Five different scenarios were assumed to obtain the soil hydraulic parameters. The HYDRUS-1D model has a database for soil hydraulic parameter estimation (uncalibrated parameters); these data were originally obtained from Rawls et al. [23] for the BC model and from ROSETTA software [47] for the VG model. Thus, the soil hydraulic parameters for the first four scenarios (S1, S2, S3, and S4) were selected from the HYDRUS-1D database of soil hydraulic properties based on the hydraulic model (BC or VG) and different soil physical information levels. To compare the uncalibrated and calibrated soil hydraulic parameters, the measured profile soil water content was used to calibrate and optimize the soil hydraulic parameters by inverse modelling (scenario S5). Finally, the obtained parameters from each previous scenario (calibrated and uncalibrated parameters) were fed into the HYDRUS-1D model to determine the accuracy of each scenario parameter by comparing the observed and simulated soil water contents. The details of the VG, BC, and other HYDRUS-1D model parameters are as follows.


### 2.2.1. Water Flow Equation

The transport of soil water was numerically simulated by solving the Richards equation (48) using the model developed by the United States Salinity Laboratory (HYDRUS-1D) [49, 50]. The HYDRUS-1D model solves the following one-dimensional Richards equation for nonsteady unsaturated water flow:(1)∂θ∂t=∂∂zKh∂h∂z+1−S,where h is the soil water pressure head [L], θ is the volumetric soil water content [L3.L−3], t is time (T), z is the soil depth [L], K(h) is the unsaturated hydraulic conductivity function [L.T−1], and S is the sink term accounting for root water uptake [L3.L−3.T−1].


### 2.2.2. Soil Hydraulic Model and Properties

To solve the Richards equation, the water retention functionθ(h) and the hydraulic conductivity function K(h) were obtained from the BC and VG hydraulic models.(1) Brooks and Corey Model. The unsaturated soil hydraulic properties θ(h) and K(h) in (1), according to Brooks and Corey [12], are given by(2)θh=θr+θs−θrαh−nh<−1α,θsh≥−1α,Kh=KsSe2/n+l+2h<−1α,Ksh≥−1α,Se=θh−θrθs−θr,where Se is the effective saturation; θr is the residual moisture content, expressed as [L3.L−3]; θh,θs,Ks, and h are as previously defined; α is the inverse of the air entry point; n is the pore size distribution index; and the HYDRUS-1D database of soil properties provides a value of l=1 for all BC soil data sets. The parameters α, n, and l are considered to be empirical coefficients that affect the shape of the hydraulic functions. As mentioned above, the BC parameters were originally reported by Rawls et al. [23].(2) Van Genuchten-Mualem Model. The unsaturated soil hydraulic properties are described by the VG model [15, 51]:(3)θh=θr+θs−θr1+αhnmh<0,θsh≥0,Kh=KsSel1−1−Se1/mm2h<0,Ksh≥0,where m=1−1/n,n>1 [52], and l was estimated by [51] as an average value of approximately 0.5 for many soil types; other variables are as previously defined. The VG parameters were obtained using ROSETTA software (VG-ROSETTA) with input data as follows: soil textural class; particle size distribution; or particle size distribution and bulk density.


### 2.2.3. Evapotranspiration

Meteorological data, such as air temperature, solar radiation, wind speed, precipitation, and relative humidity, were collected from the Changzhuang measuring station, located near the experimental site (116°51′ E, 38°21′ N), for the calibration period (4th June 2009 to 12th November 2009) and the validation period (7th May 2010 to 1st October 2010). The daily reference evapotranspiration ET0 (L.T−1) was calculated using the FAO-56 Penman–Monteith equation (53). The potential evapotranspiration rate, ETPt, was calculated by multiplying the crop coefficient (Kc) by the ET0 rate and was divided into the potential evaporation rate, EPt, and potential transpiration rate, TPt, by Beer’s Law (Ritchie [54]:(4)EPt=ETPte−σLAI,TPt=ETPt−EPt,where LAI is the leaf area index with units of [L2.L−2], and σ is a plant-dependent parameter (dimensionless), for which Ritchie [54] suggested a value of 0.4 for agricultural crops with an LAI less than 2.7 cm2 cm−2. The Kc and LAI values for a cotton crop were obtained from [55].


### 2.2.4. Root Water Uptake

The Feddes model was used to estimate the water uptake function of the root system [56].(5)Sh,z=γhβzTPt,where S(h, z) is the sink; TPt is the potential transpiration rate [L.T−1]; βz is the root density distribution function (L.L−2) (the integration of the root distribution along the root zone is equal to one; there are many ways to express the βz function, and HYDRUS-1D does not limit its users in any way); and γh is the soil water stress response function (dimensionless) of Feddes et al. [56].(6)γh=0;h≥h1,h1−hh1−h2h1>h≥h2,1;h2>h≥h3,h−h4h3−h4h3>h≥h4,0;h4≥h,.The values of the threshold parameters (h1, h2, h3, and h4) in the cotton root water uptake model were set to h1 = −1 cm, h2 = −25 cm, h3max = −200 cm, h3min = −800 cm, and h4 = −8000 cm. These values were optimized by Li [57] for cotton water uptake simulation.


## 2.2.1. Water Flow Equation

The transport of soil water was numerically simulated by solving the Richards equation (48) using the model developed by the United States Salinity Laboratory (HYDRUS-1D) [49, 50]. The HYDRUS-1D model solves the following one-dimensional Richards equation for nonsteady unsaturated water flow:(1)∂θ∂t=∂∂zKh∂h∂z+1−S,where h is the soil water pressure head [L], θ is the volumetric soil water content [L3.L−3], t is time (T), z is the soil depth [L], K(h) is the unsaturated hydraulic conductivity function [L.T−1], and S is the sink term accounting for root water uptake [L3.L−3.T−1].


## 2.2.2. Soil Hydraulic Model and Properties

To solve the Richards equation, the water retention functionθ(h) and the hydraulic conductivity function K(h) were obtained from the BC and VG hydraulic models.(1) Brooks and Corey Model. The unsaturated soil hydraulic properties θ(h) and K(h) in (1), according to Brooks and Corey [12], are given by(2)θh=θr+θs−θrαh−nh<−1α,θsh≥−1α,Kh=KsSe2/n+l+2h<−1α,Ksh≥−1α,Se=θh−θrθs−θr,where Se is the effective saturation; θr is the residual moisture content, expressed as [L3.L−3]; θh,θs,Ks, and h are as previously defined; α is the inverse of the air entry point; n is the pore size distribution index; and the HYDRUS-1D database of soil properties provides a value of l=1 for all BC soil data sets. The parameters α, n, and l are considered to be empirical coefficients that affect the shape of the hydraulic functions. As mentioned above, the BC parameters were originally reported by Rawls et al. [23].(2) Van Genuchten-Mualem Model. The unsaturated soil hydraulic properties are described by the VG model [15, 51]:(3)θh=θr+θs−θr1+αhnmh<0,θsh≥0,Kh=KsSel1−1−Se1/mm2h<0,Ksh≥0,where m=1−1/n,n>1 [52], and l was estimated by [51] as an average value of approximately 0.5 for many soil types; other variables are as previously defined. The VG parameters were obtained using ROSETTA software (VG-ROSETTA) with input data as follows: soil textural class; particle size distribution; or particle size distribution and bulk density.


## 2.2.3. Evapotranspiration

Meteorological data, such as air temperature, solar radiation, wind speed, precipitation, and relative humidity, were collected from the Changzhuang measuring station, located near the experimental site (116°51′ E, 38°21′ N), for the calibration period (4th June 2009 to 12th November 2009) and the validation period (7th May 2010 to 1st October 2010). The daily reference evapotranspiration ET0 (L.T−1) was calculated using the FAO-56 Penman–Monteith equation (53). The potential evapotranspiration rate, ETPt, was calculated by multiplying the crop coefficient (Kc) by the ET0 rate and was divided into the potential evaporation rate, EPt, and potential transpiration rate, TPt, by Beer’s Law (Ritchie [54]:(4)EPt=ETPte−σLAI,TPt=ETPt−EPt,where LAI is the leaf area index with units of [L2.L−2], and σ is a plant-dependent parameter (dimensionless), for which Ritchie [54] suggested a value of 0.4 for agricultural crops with an LAI less than 2.7 cm2 cm−2. The Kc and LAI values for a cotton crop were obtained from [55].


## 2.2.4. Root Water Uptake

The Feddes model was used to estimate the water uptake function of the root system [56].(5)Sh,z=γhβzTPt,where S(h, z) is the sink; TPt is the potential transpiration rate [L.T−1]; βz is the root density distribution function (L.L−2) (the integration of the root distribution along the root zone is equal to one; there are many ways to express the βz function, and HYDRUS-1D does not limit its users in any way); and γh is the soil water stress response function (dimensionless) of Feddes et al. [56].(6)γh=0;h≥h1,h1−hh1−h2h1>h≥h2,1;h2>h≥h3,h−h4h3−h4h3>h≥h4,0;h4≥h,.The values of the threshold parameters (h1, h2, h3, and h4) in the cotton root water uptake model were set to h1 = −1 cm, h2 = −25 cm, h3max = −200 cm, h3min = −800 cm, and h4 = −8000 cm. These values were optimized by Li [57] for cotton water uptake simulation.


## 2.3. Inverse Model Based on Uniform Design Theory

The HYDRUS-1D model has a database for soil hydraulic property estimation, and the estimated hydraulic parameters of BC and VG-ROSETTA should be accompanied by uncertainty estimates that permit an assessment of the reliability of parameter predictions. UD theory [39, 58] was used to generate the uncertainty estimates of the soil hydraulic parameters. Fang [39] made use of so-called good lattice point (glp) sets to construct many tables for UD. UD tables are denoted as follows:(7)Unqs,where U denotes the UD, n is the number of runs, s is the number of factors, and q is the number of levels.A UD table design by Fang [39] with 2 < s < 30, 3 < q < 11, and different numbers of runs (n) was used. For a large number of factors (s > 10), the number of runs is n = 110, 121, 132, 143, 154, or 165. The UD table can also be obtained from the Hong Kong Baptist University Mathematical Department website (https://www.math.hkbu.edu.hk/UniformDesign/).An uncertainty level of ±25% of the initial value of each soil hydraulic parameter was analysed; hence, each parameter was set at 11 levels (q = 11) (i.e., αinitial,αinitial±0.05,αinitial±0.10,αinitial±0.15,αinitial±0.20,andαinitial±0.25), and there was a total of five soil layers (Table 1), each of which consisted of five parameters (α, n, θs, θr, and Ks). Therefore, a total of s=5×5=25 parameters needed to be estimated. A UD table of 110 runs with 25 factors and 11 levels, designed by Fang [39], was selected U1101125. The UD test tables were obtained through the conversion of the level and actual value of each parameter. Then, the 110 sets of parameters were fed into the HYDRUS-1D model. The parameters that resulted in the highest goodness-of-fit between the observed and simulated soil water content were selected.The soil water contents were measured by time-domain reflectometry (TDR) probes on a daily basis during the calibration period (4th June 2009 to 12th November 2009) and the validation period (7th May 2010 to 1st October 2010) in soil layers at depths of 50, 100, 150, 200, 250, 300, 350, 400, and 450 cm in the study area.Table2 summarizes the input and output parameters for each scenario. For scenario S1, where the listed source is [23], parameters were created for soil texture. For scenarios where the listed source is VG-ROSETTA, parameters were derived from ROSETTA [47] using the soil textural class option for scenario S2, the particle size distribution option (sand, silt, and clay fractions) for scenario S3, and the particle size distribution and bulk density option for scenario S4. For scenario S5, the initial values of the parameters α,n,θs,θr,andKs were obtained from ROSETTA using the particle size distribution and bulk density as input data; subsequently, an uncertainty level of ±25% of the initial value of each parameter was analysed; the parameters from the UD table U1101125 that resulted in the highest goodness-of-fit between the observed and simulated soil water content were selected as the output parameters.Table 2
Input and output parameters of the model options.
ScenarioHydraulic modelSourceInput∗Output parametersParameter typeS1BC[23]Textural class (USDA)α, n, θs, θrKsUncalibratedS2VGVG-ROSETTA∗∗Textural class (USDA)α, n, θs, θrKsS3VGVG-ROSETTASand, silt, and clay fractionsα, n, θs, θrKsS4VGVG-ROSETTASand, silt, and clay fractions + soil bulk densityα, n, θs, θrKsS5VG——α, n, θs, θrKsCalibrated∗ input data for soil physical characteristics were obtained from Table 1. ∗∗ For scenarios where the listed source is VG-ROSETTA, parameters were derived from ROSETTA software [47].


## 2.4. Initial and Boundary Conditions

The initial condition was set as the initial soil water content in the study domain at timet = 0. The upper boundary condition represents the atmospheric boundary condition, where water flux depends on precipitation, evaporation, and the moisture conditions in the soil. The lower boundary condition was the measured daily pressure head derived from the groundwater table depth. The position of the groundwater table was measured every day during the calibration and validation period in an observation well located inside the experimental field (Figure 1). The ETPt values, groundwater depths, and precipitation rates during the calibration and validation periods are shown in Figure 1.Figure 1
ETPt, groundwater depth, and precipitation rate in (a) 2009 and (b) 2010.
(a)(b)


## 2.5. Statistical Analysis of the Results

Two statistical indictors were used to analyse the goodness-of-fit between the observations and simulations: (1) the root mean square error (RMSE) and (2) the Nash-Sutcliffe modelling efficiency (NSE) [59], which are calculated as follows:(8)RMSE=1N∑i=1NPi−Oi2,NSE=1−∑i=1NPi−Oi2∑i=1NOi−O¯2,where Pi represents the predicted value, Oi represents the observed value, O¯ is the mean value of all observed events at time i, and N is the number of data pairs. The RMSE provides an indication of how good the agreement between the predicted and observed values is. NSE ranged between -∞ and 1, where a high value of NSE (NSE > 0) represents a good model performance.


## 3. Results and Discussion

### 3.1. Soil Hydraulic Parameters

In this study, the soil hydraulic parameters determined by different scenarios (S1, S2, S3, S4, and S5) (Table2) were used as the input parameters in the HYDRUS-1D model to assess the reliability of the estimated parameters by comparing the simulated and measured soil water contents. Table 3 summarizes the soil hydraulic parameters for each scenario. Each scenario in Table 3 was programmed into HYDRUS-1D. Figure 2 shows a comparison of the simulated and measured soil water contents at all depths in the soil profile for the five scenarios.Table 3
Soil hydraulic parameters for the test scenarios.
ScenarioSoil depth (cm)θr (cm3.cm−3)θs (cm3.cm−3)α (cm−1)n (−)Ks (cm.day−1)l (−)Parameter typeS10–600.0150.5010.04820.21116.321Uncalibrated60–2750.0150.5010.04820.21116.321275–3450.040.4710.03070.1513.61345–4450.0410.4530.06820.32262.161445–6500.020.4370.1380.5925041S20–600.06450.43870.00511.662618.260.560–2750.06450.43870.00511.662618.260.5275–3450.09010.4820.00841.520211.110.5345–4450.03870.3870.02671.448438.250.5445–6500.0530.37470.03533.1798642.980.5S30–600.07560.44880.00551.617713.260.560–2750.0660.43730.00461.669918.110.5275–3450.08990.47630.00821.513111.290.5345–4450.03620.38920.03161.412249.040.5445–6500.04330.38380.03942.0777179.220.5S40–600.07430.43120.00561.634515.380.560–2750.06410.40620.00521.667819.180.5275–3450.08590.44680.00781.51587.480.5345–4450.03760.38340.03091.455158.80.5445–6500.04480.4710.04451.8038368.170.5S50–600.066870.452760.005881.5527812.3040.5Calibrated60–2750.067310.365580.00651.7511916.3030.5275–3450.077310.536160.008971.743175.610.5345–4450.030080.440910.038631.3823555.860.5445–6500.053760.42390.037831.53323349.7620.5Figure 2
Measured and simulated soil water contents at depths of (a) 50, (b) 100, (c) 150, (d) 200, (e) 250, (f) 300, (g) 350, (h) 400, and (i) 450 cm under scenarios S1, S2, S3, S4, and S5 in 2009.
(a)(b)(c)(d)(e)(f)(g)(h)(i)In the root zone (50 cm depth), the soil water movement efficiency based on the soil hydraulic parameters from the S1 scenario (based on the BC model) was 0.45, while values of 0.36, 0.46, and 0.20 were obtained based on the soil hydraulic parameters from the S2, S3, and S4 (based on the VG-ROSETTA model) scenarios, respectively. The efficiency for scenario S3 was higher than that for scenario S2, and the input data for scenarios S2 and S3 were soil texture and particle size distribution, respectively. However, even with the more accurate input data for S4 (particle size distribution and bulk density), the efficiency was relatively low. Such a result may have occurred because the average bulk density does not represent the bulk density distribution very well. This result was similar to that previously obtained by Garg and Gupta [60], who conducted a study where the root zone was divided into four layers (each 15 cm thick), and even in that case, a small variation in the bulk density was expected; the average efficiency for the root zone was 0.48 and 0.09 when using input data similar to those used for scenarios S3 and S4, respectively. Therefore, it is recommended that the particle size distribution be used as input data for the VG-ROSETTA model when there is limited information on the bulk density distribution and/or the bulk density varies considerably across the soil profile. The soil water movement efficiency in the root zone was 0.56 and 0.66, based on the soil hydraulic parameters from scenario S5 (calibrated parameters) in 2009 and 2010, respectively, and this scenario resulted in the best average efficiencies for all soil profiles (Table 4). Poor simulations at soil depths below 50 cm were observed from the average efficiencies for scenarios S1, S2, S3, and S4, as shown in Table 4.Table 4
NSE of soil water content simulations under different scenarios.
ScenarioYearSoil depth (cm)Average NSE50100150200250300350400450S120090.45−87.23−29.90−6.06−9.340.15−0.630.82−2.76−14.95S20.36−81.77−51.29−8.86−5.430.32−5.60−0.93−10.37−18.18S30.46−97.77−63.06−9.87−5.780.23−5.68−0.82−7.16−21.05S40.20−41.23−16.29−4.51−1.60−0.41−6.57−1.20−6.37−8.66S50.560.8340.430.830.570.680.920.860.690.71S520100.660.310.13−0.37−0.170.041.000.740.660.33It can also be observed from Table5 that the RMSE values were lower (0.004–0.038 and 0.001–0.024) for the soil water content simulations carried out using scenario S5 in 2009 and 2010, respectively, in contrast to the high RMSE values (0.018–0.115) that were obtained for the simulations conducted under scenarios S1, S2, S3, and S4. In terms of estimating the soil hydraulic parameters, scenario S5 simulated the soil water content with the highest accuracy. The statistical indicators show that the soil hydraulic parameters that were inversely determined by UD theory were accurate and reasonable.Table 5
RMSE of soil water content simulations under different scenarios.
ScenarioYearSoil depth (cm)Average RMSE50100150200250300350400450S120090.0250.0860.0430.0920.0890.0620.0460.0180.0310.05S20.0270.0830.0560.1090.0700.0560.0920.0580.0540.07S30.0250.0910.0620.1150.0720.0590.0920.0560.0460.07S40.0300.0600.0320.0820.0450.0800.0980.0620.0430.06S50.0220.0040.0060.0140.0150.0380.0100.0160.0090.01S520100.0130.0030.0080.0100.0100.0240.0010.0030.0030.01


## 3.1. Soil Hydraulic Parameters

In this study, the soil hydraulic parameters determined by different scenarios (S1, S2, S3, S4, and S5) (Table2) were used as the input parameters in the HYDRUS-1D model to assess the reliability of the estimated parameters by comparing the simulated and measured soil water contents. Table 3 summarizes the soil hydraulic parameters for each scenario. Each scenario in Table 3 was programmed into HYDRUS-1D. Figure 2 shows a comparison of the simulated and measured soil water contents at all depths in the soil profile for the five scenarios.Table 3
Soil hydraulic parameters for the test scenarios.
ScenarioSoil depth (cm)θr (cm3.cm−3)θs (cm3.cm−3)α (cm−1)n (−)Ks (cm.day−1)l (−)Parameter typeS10–600.0150.5010.04820.21116.321Uncalibrated60–2750.0150.5010.04820.21116.321275–3450.040.4710.03070.1513.61345–4450.0410.4530.06820.32262.161445–6500.020.4370.1380.5925041S20–600.06450.43870.00511.662618.260.560–2750.06450.43870.00511.662618.260.5275–3450.09010.4820.00841.520211.110.5345–4450.03870.3870.02671.448438.250.5445–6500.0530.37470.03533.1798642.980.5S30–600.07560.44880.00551.617713.260.560–2750.0660.43730.00461.669918.110.5275–3450.08990.47630.00821.513111.290.5345–4450.03620.38920.03161.412249.040.5445–6500.04330.38380.03942.0777179.220.5S40–600.07430.43120.00561.634515.380.560–2750.06410.40620.00521.667819.180.5275–3450.08590.44680.00781.51587.480.5345–4450.03760.38340.03091.455158.80.5445–6500.04480.4710.04451.8038368.170.5S50–600.066870.452760.005881.5527812.3040.5Calibrated60–2750.067310.365580.00651.7511916.3030.5275–3450.077310.536160.008971.743175.610.5345–4450.030080.440910.038631.3823555.860.5445–6500.053760.42390.037831.53323349.7620.5Figure 2
Measured and simulated soil water contents at depths of (a) 50, (b) 100, (c) 150, (d) 200, (e) 250, (f) 300, (g) 350, (h) 400, and (i) 450 cm under scenarios S1, S2, S3, S4, and S5 in 2009.
(a)(b)(c)(d)(e)(f)(g)(h)(i)In the root zone (50 cm depth), the soil water movement efficiency based on the soil hydraulic parameters from the S1 scenario (based on the BC model) was 0.45, while values of 0.36, 0.46, and 0.20 were obtained based on the soil hydraulic parameters from the S2, S3, and S4 (based on the VG-ROSETTA model) scenarios, respectively. The efficiency for scenario S3 was higher than that for scenario S2, and the input data for scenarios S2 and S3 were soil texture and particle size distribution, respectively. However, even with the more accurate input data for S4 (particle size distribution and bulk density), the efficiency was relatively low. Such a result may have occurred because the average bulk density does not represent the bulk density distribution very well. This result was similar to that previously obtained by Garg and Gupta [60], who conducted a study where the root zone was divided into four layers (each 15 cm thick), and even in that case, a small variation in the bulk density was expected; the average efficiency for the root zone was 0.48 and 0.09 when using input data similar to those used for scenarios S3 and S4, respectively. Therefore, it is recommended that the particle size distribution be used as input data for the VG-ROSETTA model when there is limited information on the bulk density distribution and/or the bulk density varies considerably across the soil profile. The soil water movement efficiency in the root zone was 0.56 and 0.66, based on the soil hydraulic parameters from scenario S5 (calibrated parameters) in 2009 and 2010, respectively, and this scenario resulted in the best average efficiencies for all soil profiles (Table 4). Poor simulations at soil depths below 50 cm were observed from the average efficiencies for scenarios S1, S2, S3, and S4, as shown in Table 4.Table 4
NSE of soil water content simulations under different scenarios.
ScenarioYearSoil depth (cm)Average NSE50100150200250300350400450S120090.45−87.23−29.90−6.06−9.340.15−0.630.82−2.76−14.95S20.36−81.77−51.29−8.86−5.430.32−5.60−0.93−10.37−18.18S30.46−97.77−63.06−9.87−5.780.23−5.68−0.82−7.16−21.05S40.20−41.23−16.29−4.51−1.60−0.41−6.57−1.20−6.37−8.66S50.560.8340.430.830.570.680.920.860.690.71S520100.660.310.13−0.37−0.170.041.000.740.660.33It can also be observed from Table5 that the RMSE values were lower (0.004–0.038 and 0.001–0.024) for the soil water content simulations carried out using scenario S5 in 2009 and 2010, respectively, in contrast to the high RMSE values (0.018–0.115) that were obtained for the simulations conducted under scenarios S1, S2, S3, and S4. In terms of estimating the soil hydraulic parameters, scenario S5 simulated the soil water content with the highest accuracy. The statistical indicators show that the soil hydraulic parameters that were inversely determined by UD theory were accurate and reasonable.Table 5
RMSE of soil water content simulations under different scenarios.
ScenarioYearSoil depth (cm)Average RMSE50100150200250300350400450S120090.0250.0860.0430.0920.0890.0620.0460.0180.0310.05S20.0270.0830.0560.1090.0700.0560.0920.0580.0540.07S30.0250.0910.0620.1150.0720.0590.0920.0560.0460.07S40.0300.0600.0320.0820.0450.0800.0980.0620.0430.06S50.0220.0040.0060.0140.0150.0380.0100.0160.0090.01S520100.0130.0030.0080.0100.0100.0240.0010.0030.0030.01


## 4. Conclusion

Uniform design theory was used to develop a new approach for an indirect method often referred to as inverse modelling in this study. The new approach is a suitable and simple tool to calibrate soil hydraulic parameters. Statistical indicators show that the soil hydraulic parameters that were inversely estimated by uniform design theory were accurate and reasonable, especially in the root layer (0–50 cm). Thus, the use of calibrated soil hydraulic parameters will increase the accuracy of simulated water movement and related agro-hydrological models.When there are no available observation data for calibrating and optimizing soil hydraulic parameters, the uncalibrated parameters obtained with the Brooks and Corey (BC) model or the van Genuchten-Mualem (VG)-ROSETTA model can be used where no a priori concept is needed. Additionally, using the BC model to estimate the soil hydraulic parameters is recommended when the soil texture is available as input data, whereas the VG-ROSETTA model is recommended when the particle size distribution is available as input data. Furthermore, when there is limited knowledge of the bulk density distribution in the soil profile, using only the particle size distribution to estimate the soil hydraulic parameters with the VG-ROSETTA model is recommended.In this study, uniform design theory with 110 runs was used to generate the uncertainty estimates of the soil hydraulic parameters. However, more soil hydraulic parameter sets could be tested by increasing the number of runs to 121, 132, 143, 154, or 165. This study obtained soil hydraulic parameters at different levels of soil information. However, another option is to measure the particle size distribution and the organic matter content, then calculate the bulk density, and extrapolate the soil water content to −33 and −1500 kPa and input this information into the VG-ROSETTA model.



---

*Source: 4544446-2022-09-30.xml*