Bayesian Inversion of Soil Hydraulic Properties from Simplified Evaporation Experiments: Use of DREAM_(ZS) Algorithm

Abstract

There is an increasing interest in identifying soil hydraulic properties from simplified evaporation experiments. However, the conventional simplified evaporation method includes a deficit due to using the linear assumption and not accounting for uncertainty in parameters. A suggested alternative method is assessing the parameter uncertainties through inverse modeling. We examined the combination of a Bayesian inverse method, namely, DREAM_(ZS), and a numerical simulation model, namely, HYDRUS-1D, for parameter inversion with data in simplified evaporation experiments. The likelihood function could be conditioned only on pressure head observations (single-objective (SO)), or on both pressure head and evaporation rate observations (multi-objective (MO)), with different treatments on the top boundary condition. Three synthetic numerical experiments were generated in terms of the soil types of sand, loam and clay to verify the inverse modeling method. The MO approach performed better than the SO approach and linear assumption when the stage 1 evaporation rate was kept constant. However, the SO inversion was more robust when oscillations existed in the potential evaporation rate. Then, the SO inverse modeling was adopted to investigate two real experiments on loamy-sand soils and compared with the linear assumption. The linear assumption could be reliable for wet conditions with stage 1 evaporation but was not always useable for a relatively dry condition, such as that with stage 2 evaporation. The inverse modeling could be more successful in capturing the whole evaporation process of soils when both stage 1 and stage 2 were involved.

Key points

• A Bayesian inverse modeling framework with DREAM_(ZS) was used for simplified evaporation experiments.
• MO inversion was more sensitive to oscillations in evaporation than SO inversion.
• Inverse modeling performed better than the linear assumption when simulating stage 2 evaporation.

1. Introduction

Soil hydraulic properties (SHPs), i.e., the dependencies of soil water content and hydraulic conductivity on the matrix suction or pressure head, are prerequisite information for quantitative analyses of hydrological processes in the vadose zone [1,2,3,4]. A variety of laboratory methods were developed to determine SHPs. Most of these methods involve a hydrostatic equilibrium condition, a steady-state flow or a transient flow experiment, such as the pressure plate extractor experiment or the one-step or multi-step outflow experiments [5,6]. These experiments provide a series of measured data of water content, pressure head and/or flow rate, which can be used to evaluate SHPs by fitting empirical curves or calibrating numerical models [7].

In the last decade, the simplified evaporation experiment [8], as one of the transient flow experimental methods that are used for the SHPs of soils, became more and more popular in laboratories, mainly due to its relative ease of implementation and ability for gaining almost continuous data points. However, the conventional simplified evaporation method (SEM) relies on various assumptions in data analyses, which caused limitations when using this method. In general, at the beginning of a simplified evaporation experiment, the saturated hydraulic conductivity is significantly higher than the flux such that the hydraulic gradient may be too small to be measured within the accuracy of pressure transducers. The criteria of the hydraulic gradient have to be considered to select available data for the relationship between the hydraulic conductivity and the pressure head [9,10,11]. Linear distributions of the pressure head and flux were conventionally assumed to estimate SHPs [8] in the SEM, leading to bias errors, particularly for coarse-textured soils [12]. Parameters of the SHPs model were generally obtained from curve fitting without considering the limited measurement range, which may include misleading errors [13]. The experiment was improved in some approaches to avoid the disadvantages of the SEM, such as using the WP4 method when observing the permanent wilting point [14,15], enhancing experimental devices to extend the measurement range and precision or using complementary additional data in the wet range [16,17,18,19].

As an alternative to the linear assumption, SHP parameters can be determined from simplified evaporation experiments through inverse modeling of the transient flow in unsaturated soil samples. The Richards equation is available to identify unbiased SHPs from evaporation experiments using inverse modeling [20,21]. Dettmann et al. [19] used both the inverse method and a linear assumption to derive SHPs in a large number of organic soil samples. In their study, the linear assumption yielded a better fitness for the wet range and the inverse method reproduced better predictions when considering the full measured range. Uncertainties and bias still exist when the optimization-based inverse approach is used [22] since they are inherent consequences of parameter estimation and model prediction of hydrological processes [3,23,24,25,26]. The quantification of uncertainties is crucial in practical implications for soil water management, designing data collection and enhancing the predictive capability of vadose zone models [27]. However, parameter uncertainties in the modeling of a simplified evaporation experiment are not well represented in the literature. To our knowledge, the only attempt was made by Minasny and Field [28], who used the generalized likelihood uncertainty estimation (GLUE) algorithm.

Inverse modeling of simplified evaporation experiments involves both uncertainties in hydraulic parameters and measurements. Parameter uncertainty in hydrological models is usually investigated through stochastic inversions with the Bayesian framework, which is widely implemented in many disciplines, such as the identification of a contaminant source [29] and the optimization of sewage management [30]. The posterior probability distribution of parameters is derived from prior knowledge, for example, using the Markov chain Monte Carlo (MCMC) algorithm [31,32]. The differential evolution adaptive metropolis algorithm (DREAM) is an efficient MCMC algorithm that was proposed by Vrugt et al. [33]; it takes advantage of a self-adaptive differential evolution learning strategy and was widely used [3,29,34,35,36,37,38]. However, this algorithm has not been introduced for the SEM besides the use of GLUE. The uncertainty in measurements may be included in the objective (likelihood) function. In previous inverse evaporation modeling experiments, both a single-objective (SO) function with measured pressure head [39] and a multi-objective (MO) function with a pressure head and evaporation flux were used [19,20,28,40]. The treatments of the upper evaporation boundary in the forward model are different between the SO and MO optimization strategies, while the performances and impacts of the metrical uncertainty have not been compared in the literature.

In this study, we developed a stochastic inverse modeling approach for estimating SHP parameters and determined their uncertainties from simplified evaporation experiments. An extended DREAM algorithm, named DREAM_(ZS) [41], was combined with HYDRUS-1D to simulate the unsaturated flow with the parameter uncertainty. The two inverse strategies (SO and MO) were performed and compared with discussions on different treatments of the upper boundary and the impacts of observation errors. This inverse modeling approach was implemented for both synthetic cases and real experiments, comparing with fitting parameters that were obtained on the basis of a linear assumption in the conventional SEM. To our knowledge, this is the first examination of DREAM_(ZS) for the parameter inversion of a simplified evaporation experiment. The difference in performance between the SO and MO approaches was found, which is worthy of attention when the DREAM_(ZS) algorithm is applied for similar experiments.

2. Materials and Methods

2.1. Simplified Evaporation Method

The SEM was proposed by Schindler [8] in which an unsteady-state vertical flow in a soil sample is driven by evaporation from the top surface that begins from a fully saturated condition. Soils are filled into a sealed container with only one open face on the top for evaporation (Figure 1a). The evaporation rate can be estimated from the measured change in the weight of the sample, while the pressure head in the soil is measured with tensiometers at two points of different heights. In general, the lower tensiometer is installed at the height of z₁ = 0.25L, while the upper tensiometer is installed at the height of z₂ = 0.75L, where L is the length of the soil sample.

The traditional SEM assumed a linear distribution of the pressure head h(z) between z₁ and z₂ (Figure 1a) so that the pressure head at the midpoint h_mid is estimated as h_mid = (h₁ + h₂)/2, where h₁ and h₂ are the measured pressure heads at the lower and upper tensiometers, respectively. The volumetric water content θ(z, t) is a function of height and time. At the midpoint, the volumetric water content θ_mid is approximated as the average water content of the soil sample, which is estimated from the measured weight loss:

$${\theta }_{\mathrm{mid}}(t)=\frac{V_0-\Delta V}{V_{\mathrm{A}}}={\theta }_s-\frac{1}{L}{\displaystyle {\int }_0^{t}E(\tau )d\tau }$$

(1)

where V₀ and ΔV are the initial water volume and reduction in water volume during the evaporation period t, respectively; V_A is the total soil sample volume; and E(τ) denotes the evaporation rate at the elapsed time τ. The pair data of θ_mid and h_mid represent the water retention function.

The relationship between the hydraulic conductivity K and the pressure head h is evaluated from the SEM by assuming a linear distribution of upward flux q(z) throughout the height. This means that the hydraulic conductivity at the middle of the soil sample K(h_mid) can be determined from the Darcy law, as follows [11]:

$$K\left(h_{\mathrm{mid}}\right)=-\frac{q_{\mathrm{mid}}}{1+(\Delta h/\Delta z)}$$

(2)

where q_mid is the water flux (upward) at the middle, which is approximated by E/2 in the SEM according to the linear assumption; Δh = (h₂ − h₁) and Δz = (z₂ − z₁); and Δh/Δz yields the pressure head gradient. The values of h_mid and Δh are calculated for the measurement interval of E with the time-averaged method. The pressure head gradient should be high enough to satisfy the accuracy of pressure transducers [11].

2.2. Numerical Modeling and Inverse Method

2.2.1. Forward Modeling of Simplified Evaporation Experiments

The vertical unsaturated flow in simplified evaporation experiments can be theoretically analyzed with the one-dimensional Richards equation as follows:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K(h)\left(\frac{\partial h}{\partial z}+1\right)\right], 0 < z< L$$

(3)

We used HYDRUS-1D [42] to numerically solve the Richards equation with specific initial and boundary conditions in simplified evaporation experiments. The initial condition was

$$h(z,t)=L-z, 0<z<L \mathrm{and} t=0$$

(4)

which referred to saturated soil under hydrostatic pressure.

At the bottom, a zero-flux boundary condition was applied:

$$K\left(\frac{\partial h}{\partial z}+1\right)=0, z=0$$

(5)

The boundary condition of the top evaporation surface can be treated with different approaches. When the SO inversion is applied, the measured evaporation rate is used to achieve a known flux boundary in the forward model [39], as follows:

$$K\left(\frac{\partial h}{\partial z}+1\right)=E(t), z=L$$

(6)

When the MO inversion is applied, the evaporation rate is an output of the forward model and is estimated as the two-stage evaporation process [11,43]:

$$K\left(\frac{\partial h}{\partial z}+1\right)=E(t), z=L \mathrm{and} E(t)=E_0, \mathrm{when} {\left.h\right|}_{z=L}<h_{\mathrm{a}}$$

(7)

$$h(z,t)=h_{\mathrm{a}}, z=L \mathrm{and} E(t)={\left.K\left(\frac{\partial h}{\partial z}+1\right)\right|}_{z=L}, \mathrm{when} {\left.h\right|}_{z=L}=h_{\mathrm{a}}$$

(8)

where E₀ is the steady evaporation rate in the first stage, which is equal to the potential evaporation that depends on atmospheric limitations, while h_a is a critical value (negative) of the pressure head at the soil surface for the second evaporation stage, which is also known as the minimum allowed pressure head in HYDRUS-1D [42]. Thus, the MO inversion introduces an additional parameter in the forward model, namely, h_a, in comparison with the SO inversion. The impact of the uncertainty about this parameter will be discussed in Section 3.2.

Empirical functions are used to parameterize SHP in simplified evaporation experiments. In HYDRUS-1D, the Mualem–van Genuchten (MvG) model [44,45] was incorporated as follows:

$$S_{\mathrm{e}}=\frac{\theta -{\theta }_{\mathrm{r}}}{{\theta }_{\mathrm{s}}-{\theta }_{\mathrm{r}}}=\left\{\begin{array}{l}{\left(1+{\left|\alpha h\right|}^n\right)}^{-m},h<0\\ 1,h\ge 0\end{array}\right.$$

(9)

$$K\left(h\right)=K_{\mathrm{s}}{S}_{\mathrm{e}}^{\eta }{\left[1-{\left(1-S_{\mathrm{e}}^{1/m}\right)}^m\right]}^2$$

(10)

where S_e is the relative saturation; θ_r and θ_s are the residual and saturated volumetric water contents, respectively; α and n are fitting parameters, where m = 1−1/n; K_s is the saturated hydraulic conductivity; and η is a pore geometry parameter.

2.2.2. Bayesian Inference with DREAM_(ZS)

To estimate the SHP parameters and determine their uncertainty from inverse modeling, we employed the extended DREAM algorithm, namely, DREAM_(ZS) [41], which is an efficient MCMC algorithm. Based on Bayesian inference, the probability density function (PDF) of the inversion parameters’ posterior distribution p(u|Y_obs) is summarized as follows:

$$p\left(u|Y_{\mathrm{obs}}\right)=\frac{p\left(u\right)p\left(Y_{\mathrm{obs}}|u\right)}{p\left(Y_{\mathrm{obs}}\right)}$$

(11)

where Y_obs denotes the observation data, u are unknown parameters listed as a vector u = (u₁, …, u_d) with the dimension size of d, p(u) represents the prior knowledge of u, p(Y_obs|u) is equivalent to the likelihood function L(u|Y_obs) with respect to observation data, p(Y_obs) is the probability of a status represented by observations that are subject to the data-generating process. p(Y_obs) is usually determined as an independent variable from p(u) [46] such that in practice, an alternative formula for the relative probability can be applied [41]:

$$p\left(u|Y_{\mathrm{obs}}\right)\propto p\left(u\right)L\left(u|Y_{\mathrm{obs}}\right)$$

(12)

For the prior parameter distribution p(u), a useful assumption when lacking prior information is the uniform distribution in a given parameter range [34]. The likelihood function is subject to both the residuals and their standard deviation when a Gaussian distribution of the residuals is assumed [2]. The standard deviation is an a priori unknown parameter and could be eliminated in the likelihood function in practice by adopting the Jeffreys prior [2,47,48], which leads to

$$L\left(u|Y_{\mathrm{obs}}\right)={\left({\displaystyle \sum _{j=1}^N{\left[y_j-f_j(u)\right]}^2}\right)}^{-\frac{N}{2}}$$

(13)

where y_j represents an observed value in measurements; f_j is the simulated result for this observation; and j = 1, …, N denotes the number of measurements.

Based on the calculated likelihood, the DREAM_(ZS) algorithm is used to generate samples from the posterior PDF. In DREAM_(ZS), multiple parallel Markov chains are run simultaneously and adaptive randomized subspaces are sampled, which speeds up convergence to the target probability distribution. The change in the position of each Markov chain (noted by i) is determined by

$$p_{\mathrm{acc}}\left(u_{\mathrm{old}}^i\to u_{\mathrm{p}}^i\right)=\min \left\{1,\frac{p\left(u_{\mathrm{p}}^i\right)}{p\left(u_{\mathrm{old}}^i\right)}\right\}$$

(14)

where $p_{\mathrm{acc}}\left(u_{\mathrm{old}}^i\to u_{\mathrm{p}}^i\right)$ is the Metropolis acceptance probability [49]. It is then compared with a random value extracted from a uniform distribution between 0 and 1. The new state is adopted for the Markov chain, i.e., $u_{\mathrm{new}}^i=u_{\mathrm{p}}^i$, when $p_{\mathrm{acc}}\left(u_{\mathrm{old}}^i\to u_{\mathrm{p}}^i\right)$ is larger. If not, the state of the Markov chain is retained, i.e., $u_{\mathrm{new}}^i=u_{\mathrm{old}}^i$. Details of the DREAM_(ZS) algorithm were presented in Vrugt [41].

As proposed by Ritter et al. [50], to avoid interference and reduce the number of inversed parameters, it is considered that some SHP parameters are not necessary to be estimated. First, we used η = 0.5 according to Mualem [44]. θ_s could be directly estimated from the weight difference of the soil sample between saturated and dried status. It was also hypothesized that h_a is a known parameter when the MO approach is applied. Thus, the unknown parameters of the inverse modeling are θ_r, α, n and K_s. A stratified random procedure, i.e., the Latin hypercube sampling (LHS) method [51], was used to generate sample values from the constrained parameter ranges (

Table 1. Reference soil parameters extracted from the HYDRUS-1D catalog to generate the three synthetic experiments.

Soil Type	θS	θr	α	n	KS	l
	(cm3 cm−3)	(cm3 cm−3)	(cm−1)	(-)	(cm min−1)	(-)
Sand	0.43	0.045 (0–0.1) †	0.145 (0–0.2)	2.68 (1–7)	0.4950 (0–1)	0.5
Loam	0.43	0.078 (0–0.1)	0.036 (0–0.2)	1.56 (1–7)	0.0173 (0–1)	0.5
Clay	0.38	0.068 (0–0.1)	0.008 (0–0.2)	1.09 (1–7)	0.0033 (0–1)	0.5

^† (a–b) denote the lower (a) and upper (b) bounds of the parameter range.

). Experiences with DREAM_(ZS) suggested that 3 parallel chains can be used to appropriately explore the posterior PDF and then avoid needing too much time for the burn-in with more parallel chains [33]. The univariate scale reduction factor, namely, the Ȓ-statistic, which compares the within-chain and between chain variances of inversed parameters, was employed to assess the convergence of the Markov chains, where Ȓ < 1.2 was adopted to indicate the convergent state [52].

2.2.3. SO and MO Inverse Strategies

Simplified evaporation experiments can provide two kinds of observation data for inverse modeling with objective functions: pressure heads at two detected positions and the evaporation rate. How the observations are included in the objective (likelihood) function is dependent on the inverse strategy.

The SO optimization strategy only uses measured pressure heads in the likelihood function. In this situation, the log-likelihood, l((u |Y_obs)) = log[L(u|Y_obs)], can be rewritten from Equation (13) as follows [41]:

$$l\left(\mathrm{u}|Y_{\mathrm{obs}}\right)=-\frac{N_h}{2}\log \left\{{\displaystyle \sum _{j=1}^{N_h}{\left[h_j^{*}-h_j\left(\mathrm{u}\right)\right]}^2}\right\}$$

(15)

where h^* and h are the measured and simulated pressure heads, respectively, and N_h is the total number of available pressure head data for the inverse modeling. The measured evaporation rate is used in the SO optimization strategy as a known flux boundary condition that is described in Equation (6).

The MO optimization strategy includes both the measured pressure heads and evaporation rate in the likelihood function, and then the corresponding log-likelihood function is given by

$$l\left(u|Y_{\mathrm{obs}}\right)=-\frac{N_h+N_E}{2}\log \left\{{\displaystyle \sum _{j=1}^{N_h}{\left[h_j^{*}-h_j\left(u\right)\right]}^2}+{\displaystyle \sum _{j=1}^{N_E}{\left[E_j^{*}-E_j\left(u\right)\right]}^2}\right\}$$

(16)

where h^* and h are the estimated normalized heads, respectively, that are estimated from the observation and numerical simulation; E^* and E are the normalized evaporation rates that are estimated from observation and numerical simulation, respectively; N_E is the total number of available evaporation data for the inverse modeling; and the number of total observations is found using N = N_h + N_E. The normalized heads and evaporation rates are calculated using the mean and standard deviation values of measurements. The simulated evaporation rate is obtained from Equations (7) and (8) with the additional parameter h_a. The potential evaporation rate, E₀, is given as a constant that is estimated from the measured evaporation data but may also include uncertainty, which will be discussed in Section 3.2.

2.3. Synthetic Numerical Experiments

In this study, synthetic numerical examples of simplified evaporation experiments were generated to check the feasibility and effectiveness of the DREAM_(ZS) based inverse modeling with the SO or MO approaches. L = 5 cm was applied for each example. Three soil textures, i.e., sand, loam and clay, were considered to generate three examples using the MvG model parameters of the SHPs listed in Table 1. The initial condition and lower boundary conditions were determined to be those presented in Section 2.2.1. The upper boundary condition was treated as being that for the two-stage evaporation, namely, Equations (7) and (8), where the potential evaporation rate, E₀, was set to 0.2 cm d⁻¹. The value of h_a was given from empirical ranges: h_a = −10³ cm for sand and loam and h_a = −10⁵ cm for clay. Nevertheless, the impact of the uncertainty in h_a values is discussed in Section 3.2. The soil column was converted to a finite element grid in HYDRUS-1D with a resolution of 0.05 cm. The convergence criteria for the moisture content and pressure head were 0.001 and 1 cm, respectively, in the numerical models. The variable time discretization was adopted when solving the Richards equation. For the clay example, the pressure head decreased rapidly during evaporation such that a small record interval, namely, 10 min, was used to produce observation data for the inverse process. The investigated duration of the clay example was determined by considering a general detectable limitation of the pressured head (−10³ cm) at the upper tensiometer, which was ~1.3 days. For the sand and loam examples, the investigated duration was extended to 14 days, similar to Iden et al. [53], and a larger record interval, namely, 100 min, was used. Synthetic observation data for the three examples are shown in Figure 1b.

3. Results and Discussion

3.1. Comparison of the Results Found Using Different Methods

The uncertainty of parameters in the three synthetic examples (sand, loam and clay) were analyzed with both inverse approaches of SO and MO from the prior parameter ranges that are listed in Table 1. Note that the h_a value used in the top boundary condition of the MO approach model was not investigated in the inverse modeling, but rather was directly given as that in Section 2.3. Inverse modeling results were obtained after 60,000 iterations with the DREAM_(ZS) algorithm along three parallel chains (each had 20,000 iterations) to ensure that a stable posterior distribution of each parameter was achieved. It required approximately 6 h of computation time (Intel(R) Core(TM) i7-10700k CPU @ 3.80 GHz; 32 GB RAM; Windows 10 Pro., 64 Bit). The difference in computational time between the SO and MO approaches was not significant.

Figure 2 shows the Markov chain trace plots of the Ȓ-statistic and parameter values for the sand example. The results of the loam and clay examples are presented in the Supplementary Materials. As indicated in Figure 2a, the DREAM_(ZS) algorithm needed about 50,000 model evaluations for the SO approach to officially reach convergence (Ȓ < 1.2). However, the MO approach needed fewer model evaluations (~20,000), as shown in Figure 2b. For the loam and clay examples, the MO approach needed more evaluations to reach convergence than that of the SO approach, while the difference in number was less than 12,000 (Figures S1 and S2). As shown in Figure 2c–f, with an increase in the number of evaluations, the samples of α and n values rapidly concentrated to a small stable range, in comparison with that of θ_r and K_s. Similar characteristics were also exhibited for the loam and clay examples (Figures S1 and S2).

Figure 3 shows the marginal posterior distributions of each sampled parameter that was obtained through the SO and MO approaches for the three synthetic examples. As indicated, the effective range of the posterior distribution was significantly smaller than the prior range. The posterior distributions of the parameters obtained from the MO approach for all the three synthetic examples sufficiently enclosed the exact value. The SO approach yielded biased posterior distributions for the sand and loam examples (Figure 3a–h); however, the maximum discrepancy was less than ~10% of the exact value and could be reduced when the time interval of the measurements was smaller than 100 min. For the clay example, the SO inversion of θ_r and K_s provided quite different values from the exact values; therefore, prior knowledge may be required for a more effective inversion. This demonstrates that both of the SO and MO approaches could approximately suggest a “unique” solution of the inverse problem when the measurements were noise-free with a high enough temporal resolution.

To compare with the traditional SEM, we also carried out a calculation using the linear assumption by fitting the θ(h) and K(h) data sets. The maximum a posteriori parameter (MAP) values that were extracted from the SO and MO approaches are listed in

Table 2. Soil hydraulic parameters that were obtained with different methods for the three synthetic examples.

Soil Type	Methods	θr	α	n	KS
		(cm3 cm−3)	(cm−1)	(-)	(cm min−1)
Sand	True value	0.045	0.145	2.68	0.495
	SO MAP	0.048	0.146	2.71	0.551
	MO MAP	0.045	0.145	2.68	0.495
	MO-T MAP	0.044	0.144	2.68	0.475
	Linear assumption	0.044	0.144	2.60	0.355
Loam	True value	0.078	0.036	1.56	0.0173
	SO MAP	0.083	0.035	1.58	0.0171
	MO MAP	0.078	0.036	1.56	0.0173
	MO-T MAP	0.100	0.026	1.76	0.0117
	Linear assumption	0.086	0.035	1.58	0.0168
Clay	True value	0.068	0.008	1.09	0.0033
	SO MAP	0.008	0.009	1.07	0.0062
	MO MAP	0.068	0.008	1.09	0.0033
	MO-T MAP	0.068	0.008	1.09	0.0033
	Linear assumption	0.000	0.007	1.08	0.0042

, which were comparable to the exact parameter values, as well as the estimation from the linear assumption (also presented in Table 2). The corresponding SHP curves of θ(h) and K(h) are depicted in Figure 4 for comparison. As shown, all the estimated SHP curves agreed well with the exact curves and are hardly distinguishable in Figure 4. An observable effect was that the data points of the linear assumption did not fall on the exact curve when pF = log(−h) was smaller than 1 (h > −10 cm), which was caused by the non-linear distribution of the flow velocity along the soil column in the early evaporation period. The difference in the performances of the methods could be more clearly identified from the parameter values given in Table 2. As indicated, the MAP value of the parameters obtained from inverse modeling through the MO approach almost fully matched the exact value, which was better than the traditional SEM result. The linear assumption also caused a deviation of about 40% when estimating the saturated hydraulic conductivity when the near-saturated data were involved. In comparison, the SO approach did not always perform better than the linear assumption.

3.2. Impacts of Critical Pressure Head for Evaporation on MO Results

A special problem that arises in the MO approach is the uncertainty of the critical pressure head for evaporation, namely, h_a. It controls the switching of the upper boundary from a prescribed flux condition (stage 1) to a prescribed head condition (stage 2), as depicted in Equations (7) and (8). Specifically, the flux on the upper boundary remains at the potential evaporation rate until the pressure head at the soil surface reaches h_a; then, the evaporation decreases but the pressure head remains at h_a. In general, h_a is empirically specified for different soil textures (higher for coarser soils), while, in theory, it is related to the air humidity and temperature [54]:

$$h_a=\frac{RT}{Mg}\ln (H_r)$$

(17)

where M is the molecular weight of water, g is the gravitational acceleration, R is the gas constant, T is the air temperature and H_r is the air humidity. An alternative way of determining h_a is to find the pressure head corresponding to the water content of θ_r + 0.005 [42].

We did not include h_a in the inversion parameters but we checked the sensitivity of the MO approach using this parameter. The inverse modeling in Section 3.1 was based on using the exact h_a value in the synthetic example. In the sensitivity analysis, we carried out additional MO inverse modeling by using a larger h_a value for each example, which was denoted as the MO-T scenario. The h_a value in MO-T was ten times the exact h_a value for the examples, i.e., −10⁴ cm in the sand and loam examples and −10⁶ cm in the clay example. This change leads to new MAP values of the parameters, which are also listed in Table 2, and new SHP curves, which are shown in Figure 4. As indicated, the changes in the results of the sand and clay examples were not significant. The clay example only experienced the stage 1 evaporation such that the inverse modeling was not affected by the h_a value. In contrast, the MO-T of the loam example yielded larger errors in the MAP value of parameters. For example, the θ_r value that was estimated from MO-T reached the upper bound of the prior range. Thus, the MO inversion approach may be influenced by the additional unknown parameter, namely, h_a, especially when stage 2 evaporation was involved. In contrast, the SO approach was not affected by this parameter because it used the actual evaporation rate as the known boundary flux in the forward modeling rather than including the evaporation rate in the likelihood function.

3.3. Sensitivity to Stochastic Oscillation in Evaporation and Observation Errors

Results in Section 3.1 and Section 3.2 were obtained when a steady-state value of the potential evaporation E₀ was applied and no errors were introduced to establish synthetic examples. In real experiments, the evaporation may be unstable due to changes in atmospheric conditions. Observation errors may also exist in the measured weight and tension of the soil sample caused by the limited accuracy of devices. In this section, we discuss the impacts of the stochastic oscillation in the evaporation rate and observation errors on the inversion results.

At first, we introduced a stochastic oscillation on the potential evaporation in the synthetic experiment of the loam example, which was synthetically determined from a stochastic Gauss noise [55]:

$$E_{\mathrm{p}}=E_0(1+\xi n_{0,1})$$

(18)

where E_p denotes the oscillated potential evaporation rate around the average value (E₀), n_0,1 is a number drawn from the standard Gaussian distribution and ξ represents the relative level of oscillation. The E₀ value, 0.2 cm d⁻¹, was equal to that in the previous synthetic numerical experiment and used to estimate E_p from Equation (18) with ξ = 0.5 for a step-by-step changing interval of 12 h. E_p is a replica of E₀ in Equation (7) in the forward modeling of the synthetic loam example. It produced a fluctuated evaporation E in stage 1 and decaying evaporation in stage 2, as shown in Figure 5a. This oscillation evaporation also led to a disturbance in the pressure heads at the lower (h₁) and upper (h₂) tensiometers in the stage 1 period, as shown in Figure 5b. However, the oscillation behavior was significantly decreased in the pressure head and did not transfer to observations in stage 2 (Figure 5c). The “observation” data in Figure 5 were then introduced for inverse modeling, while the average E_p value was applied for the MO approach, as is normally implemented. The box diagrams in Figure 6 exhibit the uncertainty of the estimated parameters around the exact value (represented by 1 in the relative form) in different scenarios. Compared to that of the no oscillation scenario (ξ = 0, Figure 6a), the ξ = 0.5 oscillation significantly increased the interquartile range of parameters that were obtained through the MO approach with a bias in the MAP value related to the exact value (Figure 6b). Inversion results obtained through the SO approach were not sensitive to the oscillation in evaporation rate. This was because only the pressure head data (with limited oscillation) were included in the likelihood function of the SO approach.

Second, we tested the effect of stochastic observation errors by using

$$Y_{\mathrm{obs}}=Y_{\mathrm{mod}}(1+\zeta n_{0,1})$$

(19)

where Y_obs is the synthetic observation data used in Equation (11), Y_mod is the simulation data in the forward modeling with or without the oscillation in evaporation and ζ represents the relative level of the measured errors. For the loam example, we used ζ = 0.2 to generate an error-including data set for inverse modeling.

As shown in Figure 6c,d, the uncertainties of the parameters were significantly enhanced by observation errors for both the SO and MO approaches. The posterior distribution of the residual water content θ_r became close to the upper bound of the prior range for the ζ = 0.2 cases. When the oscillation in evaporation did not exist, the SO approach was slightly more sensitive to observation errors than the MO approach, as indicated by the wider interquartile range of parameters (Figure 6c). Otherwise, the SO approach performed better than the MO approach by yielding a relatively smaller bias in the MAP value of the parameters (Figure 6d). The treatment of the upper boundary in the MO inversion introduced additional noises when the oscillation in evaporation existed. Among different parameters, the saturated hydraulic conductivity K_s was mostly sensitive to observation errors, while the n value was relatively insensitive to observation errors and the oscillation in evaporation. The GLUE inverse modeling of simplified evaporation experiments conducted by Minasny and Field [28] also indicated high uncertainties in K_s but this was not attributed to the evaporation oscillation and observation errors.

4. Application to Real Experiments

In the study of the synthetic numerical experiments, it was demonstrated that the Bayesian inversion of SHP with the DREAM_(ZS) algorithm could provide high-accuracy parameter values when the oscillations in the evaporation and observation errors were ignorable. In comparison to the MO approach influenced by the uncertainty of the h_a value, the SO approach was more reliable in a condition with oscillating evaporation values. In this section, we further discuss the SO approach for real experiments in comparison with the results of the linear assumption.

4.1. Soil Samples and HYPROP Experiments

Simplified evaporation experiments were carried out on two samples of loamy-sand soil, denoted as A and B. Sample A was an undisturbed soil core collected from the Ordos Plateau, China, with a saturated water content of 0.42. Sample B was filled with disturbed soils collected from the Qinghai-Tibet Plateau, China, with a saturated water content of 0.50.

The equipment, named the HYPROP system (UMS GmbH, Munich, Germany), was used to implement the simplified evaporation experiment. Tensiometers were installed at 1.25 cm and 3.75 cm underneath the soil’s surface, which were the same as the depths used in the synthetic numerical experiments. The pressure head and the weight loss of the soil samples were recorded automatically with 10 min intervals. To avoid data overrun in the inversion procedure, we extracted evaporation data of the 10 min interval and pressure head data of the 100 min interval for further analyses.

The duration of each experiment was close to 8 days. As shown in Figure 7, the evaporation of both samples experienced significant fluctuations, similar to stochastic oscillations, but only sample A experienced stage 2 evaporation. The average evaporation rates of stage 1 were approximately 0.00024 cm min⁻¹ and 0.00015 cm min⁻¹, respectively, for samples A and B.

4.2. Inversion Results

Inverse modeling of the real experiments with the SO approach also involved 20,000 model evaluations along three parallel Markov chains, where the DREAM_(ZS) algorithm needed approximately 15,000 evaluations to officially converge to a stationary distribution. We used the convergent chain states to represent the posterior distribution of the parameters.

Figure 8 shows the marginal posterior distribution of the inversion parameters. As indicated, the parameter values generally had a normal or log-normal posterior distribution. The posterior distribution of K_s for sample B seemed to be cut by the upper bound of the prior range, suggesting that more prior knowledge regarding the hydraulic conductivity may be required. The MAP value and standard deviation of the posterior parameters are listed in

Table 3. Soil hydraulic parameters of samples A and B that were obtained when using the SO inversion approach and the linear assumption.

Soils	Methods		θr	α	n	KS
			(cm3 cm−3)	(cm−1)	(-)	(cm min−1)
Sample A	Prior range		0–0.10	0–0.20	1–7	0–1
	Inversion	MAP	0.061	0.052	1.750	0.170
		S.D. †	0.003	0.002	0.021	0.021
	Linear assumption		0.086	0.017	3.406	0.228
Sample B	Prior range		0–0.15	0–0.20	1–7	0–1
	Inversion	MAP	0.122	0.014	2.020	0.770
		S.D.	0.001	0.0003	0.003	0.181
	Linear assumption		0.100	0.013	1.972	0.121

^† S.D. denotes the standard deviation of the posterior parameters.

. Cross-correlation may have existed between different inversion parameters and caused equifinality for different parameter sets. For example, θ_r and n of sample A were positively correlated with each other at a correlation coefficient that was close to 0.9. The inversion values of α and n for sample A had a negative correlation, as indicated by the correlation coefficient of −0.61. For sample B, the cross-correlation coefficients between different parameters were generally small, indicating a weaker equifinality for different parameter sets.

Estimated soil hydraulic properties were also obtained through the traditional SEM with the linear assumption. This analysis was performed by using the HYPROP-Fit software [56] and the results are also listed in Table 3. The corresponding soil water characteristic curves that were directly obtained from this method and extracted from the MAP value of inversion parameters are shown in Figure 9 for comparison. The SEM result of the θ(h) curve for sample A was quite different from the inversion result due to a strong decrease in the pressure head in the early wet period. This was also indicated by the significant difference in the α and n values of sample A that were obtained from different methods (Table 3). The θ(h) points of SEM for sample B fell close to the inversion curves. A few of the K(h) data were selected by the HYPROP-Fit, especially for sample B, and most of the data points fell below the inversion curves. The linear assumption yielded a K_s value (0.121 cm min⁻¹) for sample B that was significantly smaller than the MAP value (0.770 cm min⁻¹) of the K_s extracted from the SO inversion results.

To check the performances of the different methods, the estimated parameters were introduced into the forward modeling to reproduce changes in the pressure head at the lower (h₁) and upper (h₂) observation points. The results are shown in Figure 10, where the differences between the linear assumption and the inversion were significant for sample A and relatively insignificant for sample B. The linear assumption yielded simulated pressure heads for sample A that agreed well with observations when the evaporation time was less than about 4 days. However, it significantly underestimated the decay rate of the stage 2 evaporation. The observation data points of sample A fell into the range determined from the 95% total model uncertainties of the SO inversion, but the MAP values of the parameters underestimated the pressure head when the evaporation time was less than about 4 days. In comparison with the traditional SEM based on the linear assumption, the inverse method could be more successful in capturing the whole evaporation process of soils when both stage 1 and stage 2 evaporation were involved. The linear assumption could be more reliable for the wet condition in the stage 1 evaporation but was not always useable for a relatively dry condition, such as that in the stage 2 evaporation. Such behavior was also reported by Dettmann et al. and Iden et al. [19,53], who proposed that the fitted parameters from the linear assumption need to be evaluated when they are used in the full measured pressure head range. Our study suggested that the SO inversion was an efficient alternative method for assessing the parameters of soils from the simplified evaporation experiment.

The assumption of a uniform prior distribution (noninformative) for the inversed parameters was implied in this study, which seems to have not led to severely unreasonable estimations for the synthetic examples but could cause uncertainties in the real experiments. If prior information about the soil hydraulic parameters and their correlation are available, the parameter identifiability may be significantly improved [2]. In addition, an assumption of the used likelihood function, namely, Equation (13), is that the residuals of the model follow a Gaussian distribution. This was approximately satisfied in the synthetic examples and the real experiment of sample B but deviated in sample A. Similar non-Gaussian distributions in residuals were also reported [3]. This deviation may influence the posterior distribution of inversed parameters [57], which is worthy of further investigation using other likelihood functions.

5. Conclusions

This study implemented and examined a Bayesian inference framework for the inverse modeling of simplified evaporation experiments by combining the forward model based on HYDRUS-1D and the MCMC inversion algorithm based on DREAM_(ZS). This method has the potential to overcome disadvantages in the traditional SEM that lie on the linear assumption. There are two ways of dealing with the evaporation boundary and likelihood function. The SO inversion treats the top surface as a known flux boundary and conditions the likelihood function solely with measured pressure heads. The MO inversion treats the top surface as a convertible boundary with an additional parameter depending on the evaporation stages and conditions the likelihood function with both the observed pressure heads and evaporation rates. To examine the performance of these two approaches, three synthetic experiments were generated with HYDRUS-1D using specified parameters of sand, loam and clay, and parameters were then identified with the inversion method and linear assumption for comparison. The SO inverse modeling was finally applied to investigate two real experiments of loamy-sand soils with oscillating evaporation for almost 8 days. The following conclusions were drawn according to the obtained results:

(1)

When the evaporation in stage 1 was steady and no observation errors were introduced, the MO inversion performed better than the SO inversion and linear assumption, which yielded the values of parameters that sufficiently fell into a small range around the accurate value.

(2)

The conducted MO inversion was less sensitive to observation errors; however, the conducted SO inversion was more robust when the oscillation existed in the potential evaporation rate.

(3)

For the experiment only including the stage 1 evaporation, the SO inversion and linear assumption yielded similar parameter values and both performed well regarding reproducing time-dependent pressure heads.

(4)

For the experiment including two evaporation stages, the soil hydraulic parameters derived from the linear assumption were quite different from that identified by the SO inversion, which misestimated the time-dependent pressure head for the stage 2 evaporation.