Sensitivity and identifiability of hydraulic and geophysical parameters from streaming potential signals in unsaturated porous media

Fluid flow in a charged porous medium generates electric potentials called streaming potential (SP). The SP signal is related to both hydraulic and electrical properties of the soil. In this work, global sensitivity analysis (GSA) and parameter estimation procedures are performed to assess the influence of hydraulic and geophysical parameters on the SP signals and to investigate the identifiability of these parameters from SP measurements. Both procedures are applied to a synthetic column experiment involving a falling head infiltration phase followed by a drainage phase.

GSA is used through variance-based sensitivity indices, calculated using sparse polynomial chaos expansion (PCE). To allow high PCE orders, we use an efficient sparse PCE algorithm which selects the best sparse PCE from a given data set using the Kashyap information criterion (KIC). Parameter identifiability is performed using two approaches: the Bayesian approach based on the Markov chain Monte Carlo (MCMC) method and the first-order approximation (FOA) approach based on the Levenberg–Marquardt algorithm. The comparison between both approaches allows us to check whether FOA can provide a reliable estimation of parameters and associated uncertainties for the highly nonlinear hydrogeophysical problem investigated.

GSA results show that in short time periods, the saturated hydraulic conductivity (K_s) and the voltage coupling coefficient at saturation (C_sat) are the most influential parameters, whereas in long time periods, the residual water content (θ_s), the Mualem–van Genuchten parameter (n) and the Archie saturation exponent (n_a) become influential, with strong interactions between them. The Mualem–van Genuchten parameter (α) has a very weak influence on the SP signals during the whole experiment.

Results of parameter estimation show that although the studied problem is highly nonlinear, when several SP data collected at different altitudes inside the column are used to calibrate the model, all hydraulic $(K_{\text{s}},{\mathit{\theta }}_{\text{s}},\mathit{\alpha },n)$ and geophysical parameters (n_a,C_sat) can be reasonably estimated from the SP measurements. Further, in this case, the FOA approach provides accurate estimations of both mean parameter values and uncertainty regions. Conversely, when the number of SP measurements used for the calibration is strongly reduced, the FOA approach yields accurate mean parameter values (in agreement with MCMC results) but inaccurate and even unphysical confidence intervals for parameters with large uncertainty regions.

1 Introduction

Flow through a charged porous medium can generate an electric potential (Zablocki, 1978; Ishido and Mizutani, 1981; Allègre et al., 2010; Jougnot and Linde, 2013), called streaming potential (SP). SP signals play an important role in several applications related to hydrogeology and geothermal reservoir engineering as they are useful for examining subsurface flow dynamics. During the last decade, surface SP anomalies have been widely used to estimate aquifers' hydraulic properties (Darnet et al., 2003). Interest in SP is motivated by its low cost and its high sensitivity to water flow. Either coupled or uncoupled approaches can be used for hydraulic parameter estimation from SP signals (Mboh et al., 2012). In the uncoupled approach, Darcy velocities (e.g., Jardani et al., 2007; Bolève et al., 2009) are obtained from the tomographic inversion of SP signals and are then used for the calibration of the hydrologic model. In the coupled approach, anomalies related to the tomographic inversion are avoided by inverting the full coupled hydrogeophysical model (Hinnell et al., 2010).

The SP signals have been widely studied in saturated porous media (Bogoslovsky and Ogilvy, 1973; Patella, 1997; Sailhac and Marquis, 2001; Richards et al., 2010; Bolève et al., 2009, among others). Fewer studies focused on the application of the SP signal in unsaturated flow despite the large interest in such nonlinear problems (Linde et al., 2007; Allègre et al., 2010; Mboh et al., 2012; Jougnot and Linde, 2013). Hence, in this work we are interested in the SP signals in unsaturated porous media. Our main objective is to investigate the usefulness of the SP signals for the characterization of soil parameters. For this purpose, we evaluate the impact of uncertain hydraulic and geophysical parameters on the SP signals and assess the identifiability of these parameters from the SP measurements.

The impact of soil parameters on SP signals is investigated using global sensitivity analysis (GSA). This is a useful tool for characterizing the influential parameters that contribute the most to the variability of model outputs (Saltelli et al.,1999; Sudret, 2008) and for understanding the behavior of the modeled system. GSA has been applied in several areas, for risk assessment for groundwater pollution (e.g., Volkova et al., 2008), nonreactive (Fajraoui et al., 2011) and reactive transport experiments (Fajraoui et al., 2012; Younes et al., 2016), for unsaturated flow experiments (Younes et al., 2013), natural convection in porous media (Fajraoui et al., 2017) and seawater intrusion (Rajabi et al., 2015; Riva et al., 2015). To the best of our knowledge, GSA has never been used for SP signals in unsaturated porous media. Hence, in the first part of this study, GSA is performed on a conceptual model inspired from the laboratory experiment of Mboh et al. (2012) in which SP signals are measured at different altitudes in a sandy soil column during a falling-head infiltration phase followed by a drainage phase. Four uncertain hydraulic parameters, saturated hydraulic conductivity (K_s), residual water content (θ_r), fitting Mualem–van Genuchten parameters (α,n) and two geophysical parameters (Archie's saturation exponent (n_a) and voltage coupling coefficient at saturation (C_sat)), are investigated. GSA of SP signals is performed by computing the variance-based sensitivity indices using polynomial chaos expansion (PCE). To reduce the number of PCE coefficients while maintaining high PCE orders, we use the efficient sparse PCE algorithm developed by Shao et al. (2017) which selects the best sparse PCE from a given data set using the Kashyap information criterion (KIC).

In the second part of this study, we investigate the identifiability of hydrogeophysical parameters from SP measurements. For this purpose, parameter estimation is performed using two different approaches. The first is a Bayesian approach in which model parameters are treated as random variables and characterized by their probability density functions. With this approach, the prior knowledge about the model and the observed data is merged to define the joint posterior probability distribution function of the parameters. In the sequel, Bayesian analysis is conducted using the DREAM_(ZS) software (Laloy and Vrugt, 2012; Vrugt, 2016) based on the Markov chain Monte Carlo (MCMC) method. MCMC has been successfully used in various inverse problems (e.g., Vrugt et al., 2003, 2008; Arora et al., 2012; Younes et al., 2017). The MCMC method yields an ensemble of possible parameter sets that satisfactorily fit the available data. These sets are then employed to estimate the posterior parameter distributions and hence the optimal parameter values and the associated 95 % confidence intervals (CIs) in order to quantify the parameters' uncertainty. The second inversion approach is the commonly used first-order approximation (FOA) approach based on the standard Levenberg–Marquardt algorithm. Two scenarios are considered to check whether FOA can provide a reliable estimation of parameters and associated uncertainties for the highly nonlinear hydrogeophysical problem investigated in the case of abundant data (small uncertainty regions) and in the case of scarcity of data (large uncertainty regions). In the first scenario, SP data collected from sensors at five different locations are taken into account for the calibration. In the second scenario, only the SP data from one sensor are used for model calibration.

The present study is set out as follows. Section 2 presents the hydrogeophysical model and the reference solution. Section 3 reports on the GSA results of SP signals. Then, Sect. 4 discusses results of parameter estimation with both MCMC and FOA approaches for the two investigated scenarios.

Figure 1 Illustration of the experimental device.

2 Test case description and numerical solution

2.1 Test case description

The test case considered in this work is similar to the laboratory experiment developed in Mboh et al. (2012), involving a falling-head infiltration phase followed by a drainage phase (Fig. 1). This experiment is representative of several laboratory SP experiments (Linde et al., 2007; Allègre et al., 2010; Jougnot and Linde, 2013, among others). Quartz sand is evenly packed in a plastic tube with an internal diameter of 5 cm to a height of L_s=117.5 cm. The column is initially saturated with a ponding of L_w=48 cm above the soil surface. Five sensors allowing SP measurements are installed 5, 29, 53, 77 and 101 cm from the surface, respectively. The column has a zero pressure head maintained at its bottom. At the top of the column, the boundary condition corresponds to a Dirichlet condition with a prescribed pressure head condition during the falling-head phase, followed by a Neumann condition with zero infiltration flux during the drainage phase. During the falling-head phase, the prescribed pressure head h_top=48 cm has an exponential behavior driven by the saturated conductivity $h_{\text{top}}=\left(L_{\text{s}}+L_{\text{w}}\right)e^{-\frac{K_{\text{s}}}{L_{\text{s}}}t}-L_{\text{s}}$. The falling-head phase remains until the ponding vanishes at the critical time $t_{\text{c}}=-\frac{L_{\text{s}}}{K_{\text{s}}}\ln \left(\frac{L_{\text{s}}}{L_{\text{s}}+L_{\text{w}}}\right)$.

2.2 Mathematical model

The total electrical current density j (A m⁻²) is determined from the generalized Ohm's law as follows:

$$\begin{array}{}\text{(1)}& j=-\mathit{\sigma }\mathrm{\nabla }\mathit{\phi }+j_{\text{s}},\end{array}$$

where φ (V) is the streaming potential, j_s (A m⁻²) is the streaming current density and σ (S m⁻¹) is the electrical conductivity distribution that is assumed to be isotropic.

Hence, the conservation equation $\left(\mathrm{\nabla }\cdot j=\mathrm{0}\right)$ is written as

$$\begin{array}{}\text{(2)}& \mathrm{\nabla }\cdot \left(\mathit{\sigma }\mathrm{\nabla }\mathit{\phi }\right)=\mathrm{\nabla }\cdot j_{\text{s}}.\end{array}$$

In addition, the electrical conductivity distribution can be estimated using the saturation $S_{\text{w}}=\mathit{\theta }/{\mathit{\theta }}_{\text{s}}$ as follows (Mboh et al., 2012):

$$\begin{array}{}\text{(3)}& \mathit{\sigma }={\mathit{\sigma }}_{\text{sat}}{S}_{\text{w}}^{n_{\text{a}}},\end{array}$$

where σ_sat (S m⁻¹) is the electric conductivity at saturation and n_a is the Archie saturation exponent (Archie, 1942).

The streaming current density (j_s) can be related to the Darcy velocity (q (cm min⁻¹)) by (Linde et al., 2007; Revil et al., 2007)

$$\begin{array}{}\text{(4)}& j_{\text{s}}=\left(-{\mathit{\sigma }}_{\text{sat}}{\displaystyle \frac{\mathit{\rho }g}{K_{\text{s}}}}{C}_{\text{sat}}{S}_{\text{w}}\right)q,\end{array}$$

where K_s (cm min⁻¹) is the saturated hydraulic conductivity, ρ (kg min⁻³) is the water density, g (m s⁻¹) is the gravitational acceleration and C_sat (V Pa⁻¹) is the voltage coupling coefficient at saturation.

Hence, the combination of the previous Eqs. (1)–(4) leads to the following partial differential equation governing the SP signals:

$$\begin{array}{}\text{(5)}& \mathrm{\nabla }\cdot \left(S_{\text{w}}^{n_{\text{a}}}\mathrm{\nabla }\mathit{\phi }\right)=\mathrm{\nabla }\cdot \left(-{\displaystyle \frac{\mathit{\rho }g{C}_{\text{sat}}{S}_{\text{w}}}{K_{\text{s}}}}q\right).\end{array}$$

However, the flow through an unsaturated soil column can be modeled by the one-dimensional Richard's equation:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=\left(c\left(h\right)+S_{\text{s}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\text{s}}}}\right){\displaystyle \frac{\partial H}{\partial t}}=-\mathrm{\nabla }\cdot q,\end{array}$$

where H (cm) and h (cm) are, respectively, the hydraulic and pressure head, such that $H=h-z$; z (cm) is the depth (downward positive); S_s (–) is the specific storage; θ_s (cm³ cm⁻³) and θ (cm³ cm⁻³) are the saturated and actual water contents, respectively; c(h) (cm⁻¹) is the specific moisture capacity; and K(h) (cm min⁻¹) is the hydraulic conductivity.

The water velocity (q) is given from Darcy's law:

$$\begin{array}{}\text{(7)}& q=-K\left(h\right)\mathrm{\nabla }H.\end{array}$$

In the current study, the standard models of Mualem (1976) and van Genuchten (1980) are used to relate pressure head, hydraulic conductivity and water content:

$$\begin{array}{}\text{(8)}& \begin{array}{l}{S}_{\text{e}}\left(h\right)={\displaystyle \frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\text{r}}}{{\mathit{\theta }}_{\text{s}}-{\mathit{\theta }}_{\text{r}}}}=\left\{\begin{array}{l}{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+{\left|\mathit{\alpha }h\right|}^n\right)}^m}}h<\mathrm{0}\\ \mathrm{1}h\ge \mathrm{0}\end{array}\right.\\ K\left(S_{\text{e}}\right)=K_{\text{s}}{S}_{\text{e}}^{\mathrm{1}/\mathrm{2}}{\left[\mathrm{1}-{\left(\mathrm{1}-S_{\text{e}}^{\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}}\end{array},\end{array}$$

where S_e (–) is the effective saturation, θ_r (cm³ cm⁻³) is the residual water content, K_s (cm min⁻¹) is the saturated hydraulic conductivity, $m=\mathrm{1}-\mathrm{1}/n$, and α  (cm⁻¹) and n (–) are the Mualem–van Genuchten-shaped parameters.

2.3 Numerical model

Although several numerical techniques have been developed for the solution of the multidimensional Richards equation (e.g., Fahs et al., 2009; Belfort et al., 2009; Younes et al., 2013; Deng and Wang, 2017), the standard finite volume method is used here for the spatial discretization of the one-dimensional Richard's equation (Eq. 6). The integration of this equation over the finite volume (i) between $(i-\mathrm{1}/\mathrm{2})$ and $(i+\mathrm{1}/\mathrm{2})$ gives

$$\begin{array}{}\text{(9)}& \underset{i-\mathrm{1}/\mathrm{2}}{\overset{i+\mathrm{1}/\mathrm{2}}{\int }}\left(c\left(h\right)+S_{\text{s}}{\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\text{s}}}}\right){\displaystyle \frac{\partial H}{\partial t}}\text{d}z=q_{i-\mathrm{1}/\mathrm{2}}-q_{i+\mathrm{1}/\mathrm{2}}.\end{array}$$

Using expressions of the Darcy velocity at the element interfaces $q_{i-\mathrm{1}/\mathrm{2}}=-\frac{K_{i-\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }z}(H_i-H_{i-\mathrm{1}})$ and $q_{i+\mathrm{1}/\mathrm{2}}=-\frac{K_{i+\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{\Delta }z}(H_{i+\mathrm{1}}-H_i)$ in the case of a uniform spatial discretization with a spatial step, we obtain

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }z\left(c_i+S_{\text{s}}{\displaystyle \frac{{\mathit{\theta }}_i}{{\mathit{\theta }}_{\text{s}}}}\right){\displaystyle \frac{\partial H_i}{\partial t}}=K_{i+\mathrm{1}/\mathrm{2}}\left(H_{i+\mathrm{1}}-H_i\right)\text{(11)}& {\displaystyle }& {\displaystyle }-K_{i-\mathrm{1}/\mathrm{2}}\left(H_i-H_{i-\mathrm{1}}\right).\end{array}$$

Using $\mathit{\tau }=S_{\text{w}}^{n_{\text{a}}}$ and $\mathit{\delta }=\frac{\mathit{\rho }g{C}_{\text{sat}}{S}_{\text{w}}}{K_{\text{s}}}$, the integration of Eq. (5) over the finite volume (i) yields

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\tau }}_{i+\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }z}}\left({\mathit{\phi }}_{i+\mathrm{1}}-{\mathit{\phi }}_i\right)-{\displaystyle \frac{{\mathit{\tau }}_{i-\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }z}}\left({\mathit{\phi }}_i-{\mathit{\phi }}_{i-\mathrm{1}}\right)\text{(13)}& {\displaystyle }& {\displaystyle }-{\mathit{\delta }}_{i+\mathrm{1}/\mathrm{2}}{K}_{i+\mathrm{1}/\mathrm{2}}\left(H_{i+\mathrm{1}}-H_i\right)\text{(14)}& {\displaystyle }& {\displaystyle }+{\mathit{\delta }}_{i-\mathrm{1}/\mathrm{2}}{K}_{i-\mathrm{1}/\mathrm{2}}\left(H_i-H_{i-\mathrm{1}}\right)=\mathrm{0},\end{array}$$

where the values at the interface ${\mathit{\tau }}_{i\pm \mathrm{1}/\mathrm{2}}$, ${\mathit{\delta }}_{i\pm \mathrm{1}/\mathrm{2}}$ and $K_{i\pm \mathrm{1}/\mathrm{2}}$ are calculated using the arithmetic mean between adjacent elements (for instance, ${\mathit{\tau }}_{i+\mathrm{1}/\mathrm{2}}=\left({\mathit{\tau }}_i+{\mathit{\tau }}_{i+\mathrm{1}}\right)/\mathrm{2})$.

Then, the temporal discretization of the obtained nonlinear ODE/DAE system (9–10) is performed with the method of lines (MOL) using the DASPK (Brown et al., 1994) time solver. The MOL is suitable for strongly nonlinear systems since it allows high-order temporal integration methods with formal error estimation and control (Miller et al., 1998; Younes et al., 2009; Fahs et al., 2009, 2011). In the current study, the relative and absolute local error tolerances are fixed to 10⁻⁶.

Numerical simulations are performed assuming typical MVG hydraulic parameters for the sandy soil with (according to Carsel and Parrish, 1988) K_s=0.495 cm min⁻¹, θ_s=0.43 cm³ cm⁻³, θ_r=0.045 cm³ cm⁻³, α=0.145 cm⁻¹ and n=2.68. The voltage coupling coefficient at saturation is $C_{\text{sat}}=-\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{7}}$V Pa⁻¹ and the Archie saturation exponent is n_a=1.6.

Figure 2 Reference SP signals. Solid lines represent the reference SP solution and dots represent the sets of perturbed data serving as conditioning information for model calibration.

Based on these hydraulic and geophysical parameters, a reference (mesh-independent) solution is obtained using a uniform mesh of 235 cells of 0.5 cm length. Data are generated by sampling the output SP signals every 10 min during 1800 min. Figure 2 shows that the SP signals have an almost linear behavior in the saturated falling-head phase. During the drainage phase, they have a nonlinear behavior and approach zero voltage for the dry conditions occurring toward the end of the experiment. The SP signals are independent Gaussian random noises with a standard deviation of 2.73 10⁻⁵ V. This noise level was obtained by Mboh et al. (2012) from laboratory measurements. The noised data (Fig. 2) are used as “observations” in the calibration exercise.

3 Global sensitivity analysis of SP signals

3.1 GSA method

The aim of GSA is to assess the effect of the variation of parameters on the model output (Mara and Tarantola, 2008). Such knowledge is important for determining the most influential parameters as well as their regions and periods of influence (Fajraoui et al., 2011). The sensitivity of a model to its parameters can be assessed using variance-based sensitivity indices. These indices evaluate the contribution of each parameter to the variance of the model (Sobol', 2001). Polynomial chaos theory (Wiener, 1938) has been largely used to perform variance-based sensitivity analysis of computer models (see for instance, Sudret, 2008; Blatman and Sudret, 2010; Fajraoui et al., 2012; Younes et al., 2016; Shao et al., 2017; Mara et al., 2017). It can be stated that the PCE method is a surrogate-based approach. However, we argue that this method employs ANOVA (analysis Of variance) decomposition and hence can be considered as a spectral method (such as the Fourier amplitude sensitivity test; Cukier et al., 1973; Saltelli et al., 1999). Indeed, with this method, the sensitivity indices are directly obtained from the PCE coefficients without needing to run the surrogate model.

Let us consider a mathematical model with a random response f(ξ) which depends on (d) independent random parameters $\mathit{\xi }=\left\{{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\xi }}_d\right\}$. With PCE, f(ξ) is expanded using a set of orthonormal multivariate polynomials (up to a polynomial degree p):

$$\begin{array}{}\text{(15)}& f\left(\mathit{\xi }\right)\approx \sum _{\left|\mathit{\alpha }\right|\le p}^{s_{\mathit{\alpha }}{\mathrm{\Psi }}_{\mathit{\alpha }}\left(\mathit{\xi }\right)},\end{array}$$

where $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{1}}\mathrm{\dots }{\mathit{\alpha }}_d\in R^d$ is a dth-dimensional index. S_α represents the polynomial coefficients and Ψ_α represents the generalized polynomial chaos of degree $\left|\mathit{\alpha }\right|={\sum }_{i=\mathrm{1}}^d{\mathit{\alpha }}_i$, such as Hermite, Legendre and Jacobi polynomials, for instance. In this work, Legendre polynomials are employed because uniform distributions are considered for the parameters. The noninformative uniform distributions are used here to express the absence of prior information, which makes all possible values of the parameter equally likely.

Equation (12) is similar to an ANOVA representation of the original model (Sobol' 1993), from which it is straightforward to express V[f(ξ)], the variance of f(ξ) as the sum of the partial contribution of the inputs, as follows:

$$\begin{array}{}\text{(16)}& V\left[f\left(\mathit{\xi }\right)\right]=\sum _{\mathit{\alpha }}^{s_{\mathit{\alpha }}^{\mathrm{2}}}.\end{array}$$

The first-order sensitivity index S_i and the total sensitivity index ST_i are defined by

$$\begin{array}{}\text{(17)}& S_i={\displaystyle \frac{V\left[E\left[\left.f\left(\mathit{\xi }\right)\right|{\mathit{\xi }}_i\right]\right]}{V\left[f\left(\mathit{\xi }\right)\right]}}\in \left[\mathrm{0},\mathrm{1}\right]\end{array}$$

$$\begin{array}{}\text{(18)}& {\text{ST}}_i={\displaystyle \frac{E\left[V\left[\left.f\left(\mathit{\xi }\right)\right|{\mathit{\xi }}_{-i}\right]\right]}{V\left[f\left(\mathit{\xi }\right)\right]}}\in \left[\mathrm{0},\mathrm{1}\right],\end{array}$$

where ${\mathit{\xi }}_{-i}=\mathit{\xi }{\mathit{\xi }}_i$, E(|) is the conditional expectation operator and V(|) the conditional variance. S_i measures the amount of variance of f(ξ) due to ξ_i alone, while ST_i≥S_i measures the amount of all contributions of ξ_i to the variance of f(ξ), including its cooperative nonlinear contributions with the other parameters ξ_j. The input/output relationship is said to be additive when ${\text{ST}}_i=S_i,\forall ,i=\mathrm{1},\mathrm{\dots },d$, and in this case ${\sum }_{i=\mathrm{1}}^d{S}_i=\mathrm{1}$.

In the sequel, a PCE is constructed for each SP signal at each observable time. The number of coefficients for a full PCE representation is $P=\left(d+p\right)\mathrm{!}/\left(d\mathrm{!}p\mathrm{!}\right)$. The evaluation of the PCE coefficients requires at least P simulations of the nonlinear hydrogeophysical model. Note that P increases quickly with the order of the PCE and the number of parameters. Hence, several sparse PCE representations, in which only the significant coefficients are sought, have been proposed in the literature in order to reduce the computational cost of the estimation of the Sobol indices. For instance, Blatman and Sudret (2010) developed a sparse PCE representation using an iterative forward–backward approach based on nonintrusive regression. Fajraoui et al. (2012) developed a technique whereby only the sensitive coefficients (that affect significantly model variance) are retained in the PCE. Recently, Shao et al. (2017) developed an algorithm based on Bayesian model averaging (BMA) to select the best sparse PCE from a given data set using the KIC (Kayshap, 1982). The main idea of this algorithm is to increase the degree of an initial PCE progressively and compute the KIC until a satisfactory representation of the model responses is obtained. This algorithm is used hereafter to compute the sensitivity indices of the SP signals.

Table 1 Reference values, lower and upper bounds for hydraulic and geophysical parameters.

3.2 GSA results

The SP responses are considered for uniformly distributed parameters over the large intervals shown in Table 1. These intervals include the reference values reported in Mboh et al. (2012). The sensitivity indices of the six input parameters $(K_{\text{s}},{\mathit{\theta }}_{\text{r}}\mathit{\alpha },n,n_{\text{a}},C_{\text{sat}})$ are estimated using an experimental design formed by $N={\mathrm{2}}^{\mathrm{12}}=\mathrm{4096}$ parameter sets. The order of the sparse PCE is automatically adapted for each observable time and location. For some observable times, the PCE is highly sparse; it reaches a degree of 31 but only contains 112 nonzero coefficients.

Figure 3 Time distribution of the SP variance at 5 cm (a) and 77 cm (b) depth. The shaded area under the variance curve represents the partial marginal contributions of the random input parameters; the contribution of interactions between parameters is represented by the blank region between the shaded area and the variance curve.

Figure 3 depicts the temporal distribution of the streaming potential variance, represented by the blue curve, and the relative contribution of the parameters, represented by the shaded area. This figure corresponds to the temporal ANOVA decomposition for sensor 1 (5 cm from the soil surface) and for sensor 4 (77 cm from the soil surface). Interactions between parameters are represented by the blank region between the variance curves and the shaded area. Note that because a Dirichlet boundary condition with zero SP is maintained at the outlet boundary, the variance of the SP signal is zero at the bottom and reaches its maximum value near the soil surface. Hence, the variance is higher for the first sensor, located 5 cm from the soil surface (Fig. 3a) than for sensor 4, located 77 cm from the soil surface (Fig. 3b).

Table 2 The first-order sensitivity index S_i and the total sensitivity index ST_i for the SP signal 5 (a) and 77 cm (b) below the soil surface at different times.

The SP signals at different altitudes exhibit similar behavior (Fig. 3). In the following, we comment on the results of sensor 1 (Fig. 3a). Because K_s varies between 0.1 and 2 cm min⁻¹, the saturated falling-head phase remains until the ponding vanishes at $t_{\text{c}}=-\frac{L_{\text{s}}}{K_{\text{s}}}\ln \left(\frac{L_{\text{s}}}{L_{\text{s}}+L_{\text{w}}}\right)$. Depending on the value of K_s (see Table 1), t_c varies between t₁=20 min and t₂=403 min. Thus, in Fig. 3a, we can see that during the first time period (t≤t₁), the SP signal is strongly influenced by the value of the parameter C_sat. The first-order and total sensitivity indices at t=10 min (Table 2a) confirm that only the saturated parameters K_s and C_sat are influential. C_sat is about 17 times more influential than K_s. As expected, the remaining parameters have no influence during the first period. The total variance is 0.72 mv, and there is no interaction between the two parameters K_s and C_sat since ST_i=S_i for both of them and ${\sum }_{i=\mathrm{1}}^d{S}_i=\mathrm{1}$.

During the second period $(t_{\mathrm{1}}\le t\le t_{\mathrm{2}})$, the flow is either saturated or unsaturated depending on the value of K_s. Figure 3a shows that the variance of the SP signal exhibits its maximum value around 2.4 mv, with strong influences of the parameters K_s and C_sat and weak interactions between them (small blank region between the variance curve and the shaded area). These results are confirmed by the sensitivity indices calculated at t=70 min and reported in Table 2a for sensor 1. Both first-order and total sensitivity indices indicate that K_s is the most influential parameter. The second influential parameter is C_sat, which has a total sensitivity index about 12 times less than K_s. The parameter α is irrelevant since its total sensitivity index is 109 times less than K_s and its partial variance is $V_i=S_i\times V_i=\mathrm{0.01}$ mv, which is less than the 95 % confidence interval associated with the SP measurement (± 0.055 mv). The total variance at t=70 min is calculated to be 2.17 mv, and the output–input relationship is close to being additive since ${\sum }_{i=\mathrm{1}}^d{S}_i=\mathrm{0.94}$, which means that interactions between parameters exist but are not significant.

During the third period (t≥t₂), the variance of the SP signal reduces to 0.3 mv (Fig. 3a) and significant interactions are observed between parameters (large blank region between the shaded area and the variance curve). Table 2a shows that for t=800 min, which corresponds to dry conditions, the total variance is 0.22. First-order sensitivity indices are very small, except for θ_r. The latter is highly influential since it has a significant first-order sensitivity index (S_i=0.27) and a more significant total-sensitivity index (ST_i=0.74). The parameters C_sat and K_s are irrelevant as they have very small first-order and total sensitivity indices. Further, strong interactions are observed between the parameters since the sum of the first-order indices is far from 1 $\left({\sum }_{i=\mathrm{1}}^d{S}_i=\mathrm{0.47}\right)$. The total sensitivity indices are significantly different from first-order sensitivity indices for almost all parameters. For instance, the ratio between these two indices is around 4 for α, 5 for n_a and 7 for n. The total sensitivity index of α remains small (0.065), whereas significant total sensitivity indices are obtained for n(ST_i=0.27) and n_a(ST_i=0.47), which indicates that these two parameters are influential (although their first-order sensitivity indices are small) because of the interaction between the parameters.

Figure 3b shows similar behavior for sensor 4 located 77 cm from the soil surface. The results in Table 2b indicate that the total variance observed at t=10, 70 and 800 min are around 8 times less than for sensor 1. For the first time period, the first and total sensitivity indices are identical to those observed for sensor 1 since saturated conditions occur inside the whole column and the same effect of K_s and C_sat can be observed whatever the location inside the column. For the second time period, the sensitivity indices for sensor 4 (Table 2b) are similar to those observed for sensor 1. However, the results for the third time period show an improvement of the relevance of the parameter α, with an increase of both first and total sensitivity indices. Indeed, compared to the results of sensor 1, both first-order and total sensitivity indices tripled. Moreover, the total sensitivity index for α(ST_i=0.22) becomes close to that of n(ST_i=0.24).

In summary, the GSA applied to SP signals identifies the influential parameters and their periods of influence and shows that

• the parameter C_sat is highly influential during the first time period (t≤t₁) during which no interactions are observed between parameters;

• the parameter K_s is highly influential during the second time period $(t_{\mathrm{1}}\le t\le t_{\mathrm{2}})$ during which small interactions occur between parameters;

• the parameters θ_r,n and n_a are influential during the third time period (t≥t₂) during which dry conditions occur; during this period, strong interactions take place between parameters;

• the parameter α has no influence on the SP signals during the two first periods and presents a very small influence (S_i=0.015 and ST_i=0.065) during the third period on sensor 1 (near the surface of the column);

• the relevance of the parameter α improves with the distance from the soil surface, although the total variance diminishes with respect to this distance. The influence of α becomes significant (ST_i=0.22) on sensor 4 (located 77 cm from the soil surface) during the third period.

4 Parameters' estimation

4.1 MCMC and FOA approaches

Calibration of computer models is an essential task since some parameters (like the Mualem–van Genuchten-shaped parameters α and n) cannot be directly measured. In such an exercise, the unknown model parameters are investigated by comparing the model responses to the observations. Recently, Mboh et al. (2012) showed that the inversion of SP signals can yield an accurate estimate of the saturated hydraulic conductivity K_s, the MVG fitting parameters α and n and the Archie saturation exponent (n_a). Moreover, they showed that the quality of the estimation was comparable to that obtained from the calibration of pressure heads. In their study, Mboh et al. (2012) used the FOA approach with the shuffled complex evolution optimization algorithm SCE-UA (Duan et al., 1993).

As important as the determination of the optimal parameter sets are the associated 95 % confidence intervals (CIs) to quantify the uncertainty of the estimated values. The determination of CIs is not straightforward if the observed model responses are highly nonlinear functions of model parameters (Christensen and Cooley, 1999). In the sequel, the parameter estimation is performed using two approaches: the popular FOA approach and the Bayesian approach based on the MCMC sampler. Contrarily to FOA, the MCMC method is robust since no assumptions of model linearity or differentiability are required. Furthermore, prior information available for the parameters can be included. MCMC provides not only an optimal point estimate of the parameters but also a quantification of the entire parameter space. Several MCMC strategies have been developed for Bayesian sampling of the parameter space (Gallagher and Doherty, 2007; Vrugt, 2016). In a groundwater and vadose zone modeling context, the most widely used of these strategies is the Metropolis–Hastings algorithm (Metropolis et al., 1953; Hastings, 1970). It proceeds as follows (Gelman et al., 1996).

• i.

  Choose an initial candidate ${\mathrm{x}}^{\mathrm{0}}=\left({\mathit{\xi }}^{\mathrm{0}},{\mathit{\sigma }}^{\mathrm{0}}\right)$ formed by the initial estimate of the parameter set ξ⁰ and the hyperparameter σ⁰ and a proposal distribution q that depends on the previous accepted candidate.

• ii.

  A new candidate ${\mathrm{x}}^i=\left({\mathit{\xi }}^i,{\mathit{\sigma }}^i\right)$ is generated from the current one xⁱ⁻¹ with the generator q(xⁱ|xⁱ⁻¹) associated with the transition probability p(ξⁱ|y_mes,σ).

• iii.

  Calculate p(ξⁱ|y_mes,σ) and compute the ratio $\mathit{\alpha }=\frac{p\left({\mathit{\xi }}^i\mathrm{|}{\mathrm{y}}_{\text{mes}},\mathit{\sigma }\right)q\left({\mathrm{x}}^i\left|{\mathrm{x}}^{i-\mathrm{1}}\right.\right)}{p\left({\mathit{\xi }}^{i-\mathrm{1}}\mathrm{|}{\mathrm{y}}_{\text{mes}},\mathit{\sigma }\right)q\left({\mathrm{x}}^{i-\mathrm{1}}\left|{\mathrm{x}}^i\right.\right)}$. Additionally, draw a random number $u\in \left[\mathrm{0},\mathrm{1}\right]$ from a uniform distribution.

• iv.

  If α≥u, then accept the new candidate, otherwise it is rejected.

• v.

  Resume from (ii) until the chain $\left\{{\mathrm{x}}^{\mathrm{0}},\mathrm{\dots },{\mathrm{x}}^k\right\}$ converges or a prescribed number of iterations i_max is reached.

Many improvements have been proposed in the literature to accelerate the MCMC convergence rate (e.g., Haario et al., 2006; ter Braak and Vrugt, 2008; Dostert et al., 2009, among others). Vrugt et al. (2009a, b) developed the DREAM MCMC sampler based on the differential evolution–Markov chain method of ter Braak (2006) to improve sampling efficiency. DREAM runs multiple Markov chains in parallel and uses subspace sampling and outlier chain correction to speed up MCMC convergence (Vrugt, 2016). Laloy and Vrugt (2012) developed the DREAM_(ZS) MCMC sampler, in which a candidate for each chain is drawn from an archive of past states denoted Z, which plays the role of the generator q. Interested readers are referred to Vrugt (2016) for more details about the properties and implementation of DREAM and DREAM_(ZS). In the current study, the DREAM_(ZS) software is used for the MCMC estimation of the hydrogeophysical parameters. Note that because of the large number of model evaluations required, the MCMC method remains rarely used in practical applications compared to the FOA approach. Indeed, with FOA, the CIs are estimated once by assuming that the Jacobian remains constant within the CIs. This assumption was found to be reasonably accurate in nonlinear problems by Donaldson and Scnabel (1987). However, recently, several authors stated that parameter interdependences and model nonlinearities violate this assumption (see, for instance, Vrugt and Bouten, 2002; Vurgin et al. 2007; Gallagher and Doherty, 2007; Mertens et al., 2009; Kahl et al., 2015).

In the following, both MCMC and FOA approaches are employed for the inversion of the highly nonlinear hydrogeophysical problem using SP measurements.

4.2 Parameters estimation results

Hydrogeophysical parameters are estimated using the DREAM_(ZS) MCMC sampler (Laloy and Vrugt, 2012). Independent uniform distributions are considered for model parameter priors and likelihood hyperparameters (see Table 1). The parameter posterior distribution is written as

$$\begin{array}{}\text{(19)}& p\left(\mathit{\xi }\mathrm{|}{\mathrm{y}}_{\text{mes}},\mathit{\sigma }\right)\propto {\mathit{\sigma }}^{-N}\exp \left(-{\displaystyle \frac{\text{SS}\left(\mathit{\xi }\right)}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where SS$\left(\mathit{\xi }\right)={\sum }_{k=\mathrm{1}}^N{\left(y_{\text{mes}}^{\left(k\right)}-y_{\text{mod}}^{\left(k\right)}\left(\mathit{\xi }\right)\right)}^{\mathrm{2}}$ is the sum of the squared differences between the observed $y_{\text{mes}}^{\left(k\right)}$ and modeled $y_{\text{mod}}^{\left(k\right)}$ SP signals at time t_k for N, the total number of SP observations.

The DREAM_(ZS) software computes multiple sub-chains in parallel to thoroughly explore the parameter space. Taking the last 25 % of individuals (when the chains have converged) yields multiple sets used to estimate the updated parameter distributions and therefore the optimal parameter values and their CIs. In the sequel, the DREAM_(ZS) MCMC sampler is used with three parallel chains.

We assume that the saturated water content has been initially measured with a fair degree of accuracy. However, instead of fixing its value (as in Kool et al. , 1987, van Dam et al., 1994, and Nützmann et al., 1998, among others), we assign a Gaussian distribution to θ_s to take associated uncertainty and its effect on the estimation of the rest of the parameters into account. It is assumed here that the saturated water content was accurately measured to be θ_s=0.43 cm³ cm⁻³ by weighing the saturated soil. The corresponding error measurements are independently and normally distributed with a zero mean and a standard deviation σ_θ=0.01 cm³ cm⁻³. Hence a Gaussian distribution is assigned to θ_s, with a mean value of 0.43 cm³ cm⁻³ and a 95 % CI [0.41−0.45] cm³ cm⁻³. The rest of the hydrogeophysical parameters have noninformative uniform distributions over the ranges reported in Table 1. The error (measurement) variance is also considered to be unknown and is simultaneously estimated with the physical parameters. Two scenarios are considered to check whether the FOA approach can provide a reliable estimation of parameters and associated uncertainties for the highly nonlinear hydrogeophysical problem investigated, both in the case of abundant data (small uncertainty regions) and in the case of scarcity of data (large uncertainty regions). In the first scenario, SP data collected from the sensors located at the five locations are taken into account for the calibration. In the second scenario, only the SP data from the first sensor located 5 cm from the soil surface serve as conditioning information for model calibration. Results of the MCMC sampler are compared to those of the FOA approach for both scenarios.

Figure 4 MCMC solutions in which all SP data are considered for the calibration. The diagonal plots represent the inferred posterior probability distribution of the model parameters. The off-diagonal scatterplots represent the pairwise correlations in the MCMC drawing.

Figure 5 MCMC solutions in which calibration is performed using only SP data located 5 cm from the surface. The diagonal plots represent the posterior probability distribution of the parameters. The off-diagonal scatterplots represent the pairwise correlations in the MCMC drawing.

4.2.1 Scenario 1: inversion using all SP measurements

Figure 4 shows the results obtained with MCMC when the SP data of the five sensors are used for the calibration. The “on-diagonal” plots in this figure display the posterior parameter distributions, whereas the “off-diagonal” plots represent the correlations between parameters in the MCMC sample. Figure 4 shows nearly bell-shaped posterior distributions for all parameters. A strong correlation is observed between θ_r and n_a(r=0.98).

From the obtained MCMC sample, it is straightforward to estimate the posterior 95 % confidence interval of each parameter. This as well as the mean estimate value of each parameter obtained with both MCMC and FOA approaches are reported in Table 3.

Table 3 Estimated mean values (bold), confidence intervals (CIs) and size of the posterior CIs (italic) with MCMC and FOA approaches for scenario 1.

Table 4 Estimated mean values (bold), confidence intervals (CIs) and size of the posterior CIs (italic) with MCMC and FOA approaches for scenario 2.

The results of this table show that the parameters are well estimated from the SP measurements since (i) identified mean values are very close to the reference solution, (ii) all confidence intervals include the reference solution and (iii) the confidence intervals are rather narrow. The saturated parameters K_s and C_sat are very well estimated (with CIs around 2 %) because of data collected during the falling-head phase during which only these two parameters are influential.

The posterior CI of the parameter θ_s is similar to its prior CI. The parameter α is reasonably well estimated, with a CI around 35 %. Recall that this parameter had very small first-order and total sensitivity indices for sensor 1 but had more significant sensitivity indices for the sensors away from the soil surface (see results for sensor 4 in Table 2b). The parameter θ_r is estimated with a CI around 90 % although it was highly influential for all sensors (for instance, a first-order sensitivity index of 0.27 and a total order of 0.74 for sensor 1). The parameters n and n_a had similar GSA behavior to small first-order sensitivities (0.038 and 0.094, respectively, for sensor 1) and large total sensitivities (0.266 and 0.4715, respectively, for sensor 1); however, the inversion shows that the parameter n is well estimated with a CI less than 10 %, whereas the parameter n_a is less well estimated with a CI around 35 %. These results suggest that GSA outcomes should be interpreted with caution in the context of parameter estimation since (i) a parameter which is not relevant for the model output in one sensor can be influential for another sensor and (ii) GSA does not presume the quality of the estimation since two parameters with similar sensitivity indices can have a different quality of estimation with the inversion procedure.

Further, the results of Table 3 show that FOA and MCMC approaches yield similar mean estimated values. Moreover, very good agreement is observed between FOA and MCMC uncertainty bounds. Concerning the efficiency of the two calibration methods for this scenario, the FOA approach is by far the most efficient method since it requires only 95 s of CPU time. The MCMC method was terminated after 16 000 model runs, which required 14 116 s. The convergence was reached at around 12 000 model runs. The last 4000 runs were used to estimate the statistical measures of the posterior distribution. Recall that contrarily to FOA, MCMC can include prior information available for the parameters and allows a quantification of the entire parameter space.

4.2.2 Scenario 2: inversion using only SP measurements near the surface

In this scenario, the number of measurements used for the calibration is strongly reduced. Only SP measurements from sensor 1 (located 5 cm below the soil surface) are considered.

The results of MCMC are plotted in Fig. 5. The correlation observed between θ_r and n_a decreases slightly to r=0.95. Almost bell-shaped posterior distributions are observed for all parameters except for the parameters θ_r and α.

The results obtained with MCMC and FOA approaches depicted in Table 4 show the following.

• The FOA approach yields accurate mean estimated values similar to MCMC results for all parameters.

• The MCMC and FOA mean estimated values are close to the reference solution and to the previous scenario. The maximum difference is observed for θ_r for which the mean estimated value with scenario 2 is 15 % greater than for scenario 1.

• The MCMC CIs for the parameters K_s, θ_s, n and C_sat are close to the previous scenario. The parameters θ_s and n are well estimated (CIs < 10 %) and the parameters K_s and C_sat are very well estimated (CIs ≤ 5 %).

• Due to the reduction of the number of data used for model calibration in scenario 2, the MCMC CIs for the parameters n_a, α and θ_r are much larger than in the previous scenario. Indeed, compared to scenario 1, the CI for n_a and θ_r increases by around 60 %, whereas the CI of α is 3 times larger than for scenario 1.

• The FOA method yields accurate CIs for the parameters θ_s, n, n_a and C_sat, whereas it overestimates the CIs of θ_r (by 24 %), K_s (by 100 %) and α (by 427 %). An unphysical uncertainty region (including negative values) is obtained for the parameter α.

These results show that the FOA can fail to provide realistic parameter uncertainties and can yield larger CIs than their corresponding nonlinear MCMC counterpart. Indeed, the linearization in the FOA method assumes that the Jacobian remains constant across the CIs. This assumption was fulfilled for the first scenario in which a large number of measurements ensured small uncertainty regions. However, the assumption is not fulfilled for some parameters of the current scenario because of the large uncertainty regions induced by the reduction of the number of SP measurements.

Concerning the efficiency of the calibration methods, the FOA required approximately 174 s of CPU time, and the MCMC required many more runs to reach the convergence than in the previous scenario. Indeed, the sampler was used with 50 000 runs (35 000 runs were necessary to reach the convergence).

5 Conclusions

In this work, a synthetic test case dealing with SP signals during a drainage experiment has been studied. The test case is similar to the laboratory experiment developed in Mboh et al. (2012), involving a falling-head infiltration phase followed by a drainage phase. GSA and Bayesian parameter inference have been applied to investigate (i) the influence of hydraulic and geophysical parameters on the SP signals and (ii) the identifiability of hydrogeophysical parameters using only SP measurements. The GSA was performed using variance-based sensitivity indices which allow the contribution of each parameter (alone or by interaction with other parameters) to the output variance to be measured. The sensitivity indices have been calculated using a PCE representation of the SP signals. To reduce the number of coefficients and explore PCE with high orders, we used the efficient sparse PCE algorithm developed by Shao et al. (2017), which selects the best sparse PCE from a given data set using the Kashyap information criterion (KIC).

The GSA applied to SP signals showed that the parameters C_sat and K_s are highly influential during the first period corresponding to saturated conditions. The parameters θ_r, n and n_a are influential when dry conditions occur. In such conditions, strong interactions take place between these parameters. The parameter α has a very small influence on the SP signals near the soil surface but its sensitivity increases with depth although the total variance decreases with depth.

Parameter estimation has been performed using MCMC and FOA approaches to check whether FOA can provide a reliable estimation of parameters and associated uncertainties for the highly nonlinear hydrogeophysical problem investigated. All hydraulic (K_s, θ_r, α and n) and geophysical (n_a and C_sat) parameters can be reasonably estimated in the first scenario for which the whole SP data (measured at five different locations) are used as conditioning information for the model calibration. The confrontation with GSA results shows that the latter should be interpreted with caution when used in the context of parameter estimation since (i) a parameter which is not relevant for the model output in one sensor can be influential for another sensor and (ii) GSA does not presume the quality of the estimation since two parameters with similar sensitivity indices can have a different quality of estimation with the inversion procedure (see, for instance, parameters n and n_a). Furthermore, although the studied problem is highly nonlinear, the FOA approach provides accurate estimations of both mean parameter values and CIs in the first scenario. These results are identical to those obtained with MCMC.

When the number of SP measurements used for the calibration is considerably reduced (i.e., data are scarce), the MCMC inversion provides larger uncertainty regions of the parameters. The FOA approach yields accurate mean parameter values (in agreement with MCMC results) but inaccurate and even unphysical CIs for some parameters with large uncertainty regions.