Coupled inversion of hydraulic and self‐potential data from transient outflow experiments to estimate soil petrophysical properties

Hydraulicproperties of soils could play an important role in affecting the partitioning of precipitation in the critical zone. In addition to traditional approaches, in the last two decades, many geophysical methods have been used to aid the hydrologic characterization and measurement of geological materials. In particular, the self‐potential (SP) method shows great potential in these hydrogeophysical applications. The objective of this study is to evaluate whether the addition of SP data can improve the estimation of hydraulic properties of soils in an outflow experiment. A stochastic, coupled hydrogeophysical inversion was developed, in which the governing equations were solved using the finite volume method and the parameter estimation was conducted using a Bayesian approach associated with the Markov chain Monte Carlo technique. The results show that the addition of SP data in the inversion could reduce the uncertainty related to the estimated hydraulic parameters of soils and the length of the associated 95% confidence interval can be shortened by ∼1/3. It is also shown that the electrical properties of soils at saturated and unsaturated conditions may also be estimated from the outflow experiment when SP data are available. Compared with hydraulic parameters, the accuracy of the estimated electrical properties is slightly lower. Among them, the saturated streaming potential coupling coefficient Csat has the highest accuracy and lowest uncertainty since Csat directly influences the magnitude of SP signals. The accuracy of other electrical parameters is lower than that of Csat (and hydraulic parameters), and the associated uncertainty can be one order of magnitude larger.


INTRODUCTION
The spatial variability of hydraulic properties in the subsurface could significantly affect the partitioning of precipitation in the critical zone (Takagi & Lin, 2011), an open system extending from the top of the canopy to the base of active portion of the critical zone, the vadose zone plays a paramount role in affecting the water chemistry and dynamics in the critical zone (Jin et al., 2011). In particular, hydraulic properties of the vadose zone such as soil water retention curve and unsaturated hydraulic conductivity regulate the amount of water stored in the soil and drained into fractured bedrock and streams (Hammond et al., 2019) and thus could be the dominating parameters influencing the hydraulic partitioning in the critical zones in, for example, mountainous western United States. In addition to critical zone hydrology, the hydraulic properties of unsaturated soils are also important in other disciplines such as agriculture engineering (Kirkham, 2014), geotechnical engineering (Lu, 2020), and environmental sciences (Leme & Miguel, 2018). In general, hydraulic properties of unsaturated soils can be either estimated with pedotransfer functions (Wösten et al., 2001) or measured with laboratory or field hydraulic tests (Tarantino et al., 2009). Pedotransfer functions, which are usually developed based on statistical regression analysis or theoretical modeling, predict soil hydraulic properties from some easy-to-measure properties usually available from soil surveys (Aimrun & Amin, 2009;Bouma, 1989;Jaiswal et al., 2013;Patil & Singh, 2016). Despite their simplicity, pedotransfer functions may not perform well when the chemical composition, texture, and/or fabric of the soil exceed the calibration range (Patil & Singh, 2016). In the laboratory, the unsaturated soil properties can be accurately measured with some steady-state hydraulic experiments such as the constant flow method proposed in Lu et al. (2006). However, this type of test is usually time consuming, especially for fine materials with low hydraulic conductivity (Moebius et al., 2007;Šimůnek et al., 1998). In contrast, a transient flow test could produce large amounts of hydraulic data (e.g., pore water pressure and flow rate) within a short time period (Gribb, 1996;Wayllace & Lu, 2012). Unsaturated hydraulic properties of the sample can be estimated from these transient data by either analytical calculations (Li et al., 2009) or inversion. In such an inversion, the difference between the measured and simulated transient hydraulic responses is minimized (Latorre et al., 2015;Šimůnek et al., 1998) to recover the soil hydraulic properties. This inversion-based parameter estimation can be either deterministic (Šimůnek & van Genuchten, 1997) or stochastic (Thoma et al., 2014). Compared with steadystate tests, transient tests reduce the measurement time significantly, and thus they have been increasingly adopted in practice to determine the unsaturated hydraulic properties of various soils (Bahrami & Aghamir, 2020;Elliott & Price, 2020).
Recently, geophysical measurements have been used to aid the estimation of hydraulic properties and to monitor various hydrological processes in the subsurface . These hydrogeophysical applications are based on the observed correlations between hydrological and geophysi-

Core Ideas
• Self-potential data can reduce the uncertainty of hydraulic parameter estimation. • Soil electrical properties can be estimated from transient hydraulic and SP data. • Stochastic coupled inversion is suitable for hydrogeophysical parameter estimation.
cal properties of porous media (Lesmes & Friedman, 2005), such as water content and dielectric constant (Topp et al., 1980), hydraulic conductivity and resistivity (Purvance & Andricevic, 2000), and hydraulic conductivity and imaginary conductivity (Weller et al., 2015). One of the most prominent geophysical methods in hydrological applications is the self-potential (SP) method, which measures the natural occurrence of electric fields on the ground surface (Allègre et al., 2010;Mboh et al., 2012). The dominant contribution to the SP signals in hydrological settings is the so-called streaming potential, which is generated by water flow in geological materials with charged mineral surfaces (Sill, 1983;Revil et al., 2002). In general, the measured streaming potential is influenced by the water flow rate, surface charge density, and effective electrical conductivity of the material. Due to the direct coupling between water flow and streaming potential, the SP method has been used in groundwater hydrology (Revil & Jardani, 2013), for example, to estimate the hydraulic conductivity and geometry of an aquifer (Darnet et al., 2003), to locate flow pathways and estimate the seepage velocity in dams and embankments (Bolève et al., 2009), and to monitor transpiration-induced water flows (Voytek et al., 2019) and rainwater infiltration processes (Doussan et al., 2002;Hu et al., 2020;Jougnot et al., 2015) in the vadose zone. At saturated conditions, the streaming potential of porous media has been extensively studied (Ishido, 1989), and theoretical models are available to describe or predict the streaming potential coupling coefficient (Hunter et al., 2013), a parameter quantifying the streaming potential generated in a material under a given pressure gradient. For unsaturated porous media, many models have also been developed to describe the coupling coefficient; these unsaturated models are developed either by modifying the Helmholtz-Smoluchowski equation (Darnet & Marquis, 2004;Guichet et al., 2003;Perrier & Morat, 2000;Thanh et al., 2020) or by upscaling the effective excess charge defined at the pore scale Revil et al., 2007). The upscaling of the effective excess charge can be done by the volume averaging approach  and flux averaging approach (Jougnot et al., 2012;Soldi et al., 2019). It should be noted that existing unsaturated coupling coefficient models may work for one type of soils but fail for a different soil type. For instance, the model developed in Linde et al. (2007) has been successfully applied to sand (Jougnot & Linde, 2013;Mboh et al., 2012), but it fails to describe the SP responses of other soils (Jougnot et al., 2012Zhang et al., 2017).
Recently, the SP method has been used in the laboratory to monitor the transient water flow in soil column tests (Allègre et al., 2014;Jougnot & Linde, 2013;Linde et al., 2007;Mboh et al., 2012). Due to the advent of coupled hydrogeophysical inversion schemes (Hinnell et al., 2010), it becomes possible to recover the unsaturated hydraulic properties of soils from both hydraulic and SP data (Mboh et al., 2012;Younes et al., 2018). There are two potential yet unproven benefits of incorporating SP measurement in a transient outflow experiment. First, SP signals induced by water flow could help better constrain some hydraulic parameters (e.g., saturated hydraulic conductivity; Schwärzel et al., 2006) that usually have high uncertainty if inverted from transient hydraulic data alone. Second, in addition to hydraulic properties, the electrical properties of soils in both saturated and unsaturated conditions may also be estimated (Younes et al., 2018). Note that accurately measuring the electrical properties of unsaturated soils is a challenging task because special instruments are often required to maintain a stable unsaturated condition (Merritt et al., 2016;Wu et al., 2017). Considering that many outflow experiment systems have been equipped with SP monitoring electrodes (Allègre et al., 2014;Linde et al., 2007;Mboh et al., 2012), it is thus necessary to quantitatively evaluate whether the estimation of hydraulic and electrical properties will benefit from the incorporation of SP data.
The objective of this study is to conduct such an evaluation using hydraulic and SP monitoring data during outflow experiments. These data will be inverted in the Bayesian framework with the coupled hydrogeophysical inversion approach (Hinnell et al., 2010) to estimate the hydraulic and electrical properties of soils. In this paper, we first introduce the coupled forward modeling of water flow and streaming potential generation in saturated and unsaturated soils, followed by the stochastic, coupled inversion, which adopts the adaptive Metropolis (AM) algorithm to estimate model parameters as well as their uncertainties. Synthetic outflow experiments are conducted on sand and loam samples to produce time-series cumulative flow, water pressure, and SP data. These datasets are inverted to obtain the electrical and hydraulic properties of soils, which are then compared with the true values to evaluate the benefits of including SP data. The stochastic, coupled inversion is also performed on published experimental data. Discussions and major conclusions are presented at the end of the paper.

FORWARD MODELING OF WATER FLOW AND STREAMING POTENTIAL
In this section, we describe the governing equations for onedimensional water flow and streaming potential in saturated and unsaturated soils. Constitutive relationships describing hydraulic (water content and hydraulic conductivity) and electrical properties (electrical conductivity and streaming potential coupling coefficient) of soils under different pore water pressures or degrees of saturations are also introduced. In this study, the governing equations are solved numerically with the finite volume method, and the details of the calculation are also presented in this section.

Water flow in saturated and unsaturated soils
Water flow in saturated and unsaturated soils can be described by the Richards equation (Richards, 1931). Under defined initial and boundary conditions, the transient hydraulic responses of the soil can be determined by solving the Richards equation (van Dam & Feddes, 2000). The variable to be solved can be either water content or pore water pressure (head). The Richards equation, which combines Darcy's law and mass conservation (Brunone et al, 2003;Namin & Boroomand, 2012), is a highly nonlinear partial differential equation (Caviedes-Voullième et al., 2013), expressed as where θ is the volumetric water content (m 3 m −3 ), h is the pore water pressure head (m), z is the vertical coordinate (positive upward, m), t is time (s), and K is the hydraulic conductivity (m s −1 ). Note that both θ and K are not constant but functions of h, known as the soil water retention curve θ(h) and hydraulic conductivity function K(h). In this study, we use the van Genuchten-Mualem models to describe θ(h) and K(h), which can be expressed respectively as (Caviedes-Voullième et al., 2013;Mualem, 1976;van Genuchten, 1980) and where θ s and θ r denote the saturated and residual water content (m 3 m −3 ), respectively, K s is the saturated hydraulic conductivity (m s −1 ), n is a unitless parameter characterizing the shape of the curve θ(h) and is mainly influenced by the pore size distribution of the material, α is a fitting parameter interpreted as the inverse of the air-entry pressure (m −1 ), and S e is the effective saturation expressed as e = θ − θ r θ s − θ r (4)

Streaming potential
In a streaming potential problem, the total electrical current density j in a porous medium can be expressed as the sum of the conductive current density j c and streaming current density j s (Sill, 1983): = + The conductive current density j c is related to the electrical potential ϕ by Ohm's law: where σ is the effective electrical conductivity of the medium (assumed to be isotropic, S m −1 ). Applying the continuity condition (i.e., ∇ ⋅ = 0) to Equation 5 results in the governing equation for the streaming potential: The effective electrical conductivity of porous geological media is influenced by many factors, including soil texture, water content, pore water chemistry, and mineral surface properties (Friedman, 2005). In this study, the empirical Archie's law is used to describe σ (Friedman, 2005): where σ w is the electrical conductivity of water (S m −1 ), F = θ s -m is the formation factor (m being the porosity exponent or cementation exponent), S = θ/θ s is the degree of saturation, n a is the Archie saturation exponent (Archie, 1942), and σ s represents the surface conductivity (S m −1 ; Revil & Glover, 1998). Equation 8 assumes that the surface conduction is in parallel with the conduction contributed from bulk water. Other more sophisticated electrical conductivity models (Bussian, 1983) may also be used here.
The streaming current density j s can be explained by the electrokinetic theory (Ishido & Mizutani, 1981). In geologi-cal materials, the surface of minerals is usually charged, and ions in the pore water can accumulate near the solid-liquid interface in response to the charged mineral surface, forming the electrical double layer (EDL; Revil & Cerepi, 2004;Titov et al., 2005). The movement of water in the pore space will drag a portion of the excess charges in EDL, resulting in the streaming current (Revil & Jardani, 2013). In principle, the streaming current density in a porous medium is related to the amount of excess charge moving with water (i.e., the excess charge outside the shear plane; see Jougnot et al., 2020) and the velocity of the pore water (Leroy & Revil, 2004).
Considering a small length Δl of the material along the flow direction, the streaming current-induced electrical potential difference across this length Δϕ can be determined using Ohm's law: where j s is the magnitude of the streaming current in the flow direction. The streaming potential coupling coefficient C is defined as the ratio of where ρ w is the water density (kg m −3 ) and g is the gravity acceleration (m s −2 ). Consider Darcy's law for the length Δl and then the Darcy velocity u in the flow direction can be expressed as Inserting Equations 10 and 11 into Equation 9 yields the relationship between the streaming current density and Darcy velocity, which is expressed as Revil et al., 2007;Younes et al., 2018) If the Darcy velocity is treated as a vector (i.e., u), the vector form of the streaming current density can be expressed as Inserting Equation 13 into Equation 7 gives the governing equation for the streaming potential incorporating C: It is noted that the parameters K and σ are dependent on the degrees of saturation (e.g., Equation 3 for K and Equation 8 for σ). The coupling coefficient C is also a function of saturation, and it can be linked to the coupling coefficient at saturation C sat , for example, by where σ sat is the effective electrical conductivity of the soil at saturation (i.e., Equation 8 with S = 1). Note that although Equation 15 has been successfully applied to sand (Jougnot & Linde, 2013), it may not work well for other types of soils (Zhang et al., 2017). Other constitutive models for C may be used (Perrier & Morat, 2000) to model the streaming potential in unsaturated soils when necessary.

Finite volume method for the coupled forward modeling
For coupled water flow and streaming potential modeling, the governing equations include Equations 1 and 14 and the relevant constitutive models include Equations 2, 3, 8, and 15. Due to the highly nonlinear dependence of θ, K, σ, and C on h, analytical solutions of Equations 1 and 14 may only be available under some specific boundary conditions and simplified constitutive models (Namin & Boroomand, 2012). In general conditions, numerical methods are required to solve Equations 1 and 14 to determine the spatial and temporal distributions of h and ϕ in the soil. In this study, we use the finite volume method (Caviedes-Voullième et al., 2013;Eymard et al., 1999) to solve the governing equations, and many similar applications can be found in the literature (Manzini & Ferraris, 2004;Namin & Boroomand, 2012;Pei et al., 2006;Younes et al., 2018).
In our finite volume modeling, the sample domain is evenly divided into N − 1 layers (with a thickness of ∆z) by N notes. Define the pressure head and water content at node i (i = 1, . . . , N) as θ i and h i , and then the discretized Equation 1 at time t reads where K i+1/2 (or K i − 1/2 ) represents the interlayer hydraulic conductivity between layer i and i + 1 (or i − 1). The interlayer conductivity can be calculated as the arithmetic mean, geometric mean, harmonic mean, or upstream mean of adjacent two layers (Baker, 2006;Caviedes-Voullième et al., 2013;Pei et al., 2006;Warrick, 1991). In general, when the hydraulic conductivity does not change significantly between adjacent layers, the simple arithmetic mean should be sufficient; otherwise, the other methods may be used to ensure numerical accuracy. In this study, the arithmetic mean was adopted [i.e., Similarly, define the electrical conductivity and electrical potential at node i as σ i and ϕ and parameter δ = (ρ w gσC)/K at node i as δ i = (ρ w gσ i C i )/K i . The discretized Equation 14 can then be written as where σ i + 1/2 (or σ i − 1/2 ) is the interlayer electrical conductivity between layers i and i + 1 (or i − 1) and δ i + 1/2 (or δ i − 1/2 ) is the average δ of two adjacent layers i and i + 1 (or i − 1). In this study, the arithmetic mean is used to calculate σ i The time domain to be solved can be divided into M − 1 steps with increment Δt j where j = 1, . . . , M. The time step Δt j may vary, and in this study, we use the implicit scheme in Caviedes-Voullième et al. (2013) to ensure that the solutions are accurate at large time steps. The discretized Equation 1 in the time domain can be expressed as (Caviedes-Voullième et al., 2013;Kumar, 1996) where the index j indicates the associated variable at the time node j and the index j + 1/2 indicates the average soil property during the time step Δt j . In this study, the arithmetic mean is used to calculate the average K during a time increment. For Equation 14, temporal discretization is not needed, and Equation 17 will be solved after Equation 18 is solved at each time node. After the water content of the soil sample is calculated, the related cumulative outflow increment at the time node j (i.e., ΔQ j ) can then be determined by considering the average water content change of the sample. The cumulative outflow at time node j (i.e., Q j ) thus can be expressed as where A is the cross section area of the sample.
In the finite volume modeling of water flow, the following boundary conditions (Dirichlet and Neumann types) are applied (Allègre et al., 2010): Vadose Zone Journal where h b (t) and u b (t) are the pressure head and Darcy velocity at the boundaries. In streaming potential modeling, Dirichlet and Neumann type boundary conditions are also applied, expressed respectively as and where ϕ b and j b are the electrical potential and external current density at the boundary.

STOCHASTIC COUPLED INVERSION
The hydraulic and geophysical properties of soils at saturated and unsaturated conditions can be estimated from the transient hydrogeophysical data using Bayesian inference, which uses probability distributions to describe model parameters m. Assume the prior information on the model parameters can be represented by the prior probability distributions P 0 (m). In Bayesian inference, the prior distributions P 0 (m) are updated as more information (e.g., measurement data d) becomes available, yielding the posterior probability distributions π(m|d) for the model parameters. From π(m|d), we can obtain not only a maximum likelihood model but also a quantification of the uncertainties of the model. Using Bayesian methods for parameter estimation has received increased popularity in many earth and environmental sciences such as atmosphere (Smith et al., 2009;Tamminen, 2004), geophysics (Grana et al., 2017;Ray & Myer, 2019), hydrology (Freni & Mannina, 2010;Tang et al., 2016), and environmental sciences (Ahmadi et al., 2015;Liu et al., 2021). Due to the analytically intractable nature of many forward modeling problems, the implementation of the Bayesian method for parameter estimations is usually aided by the Markov chain Monte Carlo (MCMC) techniques, which use random walk approaches to generate samples that follow the posterior distributions of the model parameters (Vrugt et al., 2003). Over the years, many MCMC algorithms have been developed to sample the model space, including the Metropolis-Hastings algorithm (Hastings, 1970;Metropolis et al, 1953) and its variants such as the delayed rejection adaptive algorithm (Tierney & Mira, 1999). In this study, we adopted the AM algorithm (Haario et al., 2001) to sample the model parameters that characterize the hydraulic and electri-cal properties of soils under saturated and unsaturated conditions.

Adaptive Metropolis algorithm
The AM algorithm is based on the conventional Metropolis algorithm with symmetric Gaussian proposal distributions, and during the MCMC sampling process, the sizes and orientations of the proposal distributions vary (Haario et al., 2001).
The AM algorithm has the advantages of keeping detailed balance and ergodicity and showing great efficiency on complex and highly nonlinear target distributions (Saksman & Vihola, 2010;Tamminen, 2004). In practice, the AM algorithm can be realized with the following steps (Tamminen, 2004). Assume we have already sampled k model vectors m 0 , . . . , m k − 1 . To get the next model vector m k , a candidate vector Z is sampled from the Gaussian proposal distribution with the mean value at the current point m k − 1 and with the covariance matrix C k . The covariance matrix can be expressed as where K k = cov (m 0 , . . . , m k −1 ), s n = (2.4) 2 /n is a scaling parameter (n being the dimension of the vector Z), ε is a small value (e.g., 10 −10 ) to prevent C k from being singular, I is an n × n identity matrix, and k 0 defines the burn-in period, during which, the covariance matrix C k is not updated. The samples in the burn-in period will be discarded in the calculation of the posterior distributions such that the impact of initial point (m 0 ) can be minimized (Tamminen, 2004).
In the next step, the sampled vector Z is either accepted (i.e., m k = Z) or rejected (thus, m k = m k − 1 ). The probability of accepting Z is In the AM algorithm, the ratio π( )∕π( −1 ) can be expressed as (Tamminen, 2004) where the likelihood P(d|Z) [or P(d|m k − 1 )] is a measure of the degree of fit between observed data d and data predicted from forward modeling with parameters Z (or m k − 1 ) (e.g., see Mosegaard & Tarantola, 1995) and the prior distribution P 0 (Z) [or P 0 (m k − 1 )] contains our existing knowledge of the model parameters. Thus, the computation of the acceptance F I G U R E 1 Flowchart of the coupled stochastic inversion using Markov chain Monte Carlo (MCMC) method. Q is the cumulative outflow; h is the pore water pressure head; θ is the volumetric water content; ∆V is the measured self-potential (SP); m k and m k − 1 are model vectors at step k and k − 1, respectively; Z is a candidate vector; P(d|Z) [or P(d|m k − 1 )] is the likelihood measuring the degree of fit between observation d and Z (or m k − 1 ); and β(m k − 1 , Z) is the probability of accepting Z ratio involves mainly the evaluation of the likelihood. The above steps are iterated until k reaches a predefined number.

MCMC inversion of transient hydrogeophysical data
In this study, the AM algorithm is used to obtain sample realizations of the model vector m. The associated parameters include θ r , α, n, log(K s ), m, n a , and C sat × 10 −7 . Using log(K s ) is to ensure the positivity of K s in the inversion and multiplying C sat by 10 −7 will scale up the coupling coefficient to the same order of magnitude as the other parameters. Other parameters in the constitutive models such as θ s and σ w are relatively easy to determine in practice, and thus they are assumed known and will not be estimated in our inversion. The transient hydraulic data include the pore water pressure head h, cumulative outflow Q, and SP ϕ.
In this study, the hydraulic data and electrical data are inverted in a coupled way, which relies on the direct coupling between the streaming current j s and Darcy velocity u (Equation 13) in the forward modeling and inversion process (Hinnell et al., 2010). The workflow of coupled hydrogeophysical inversions has been detailed in Ferré et al. (2009) and is also shown in Figure 1. In such a coupled inversion, an initial model is proposed, used to simulate transient responses based on coupled hydrologic and geophysical simulations (Equations 1 and 14), and then updated based on misfit between the simulated and observed transient data. In addition to the streaming current, the electrical conductivity of the soil is also dependent on the water content (e.g., Equation 8). In the inversion, both geophysical responses (ϕ) and hydraulic responses (Q and h) are used in quantifying the misfit between predicted and measured observations. Comparing to sequential inversion (Kang et al., 2020), the coupled inversion does not involve an intermediate geophysical inversion step before conducting the hydraulic inversion (or estimation), and thus it has the potential to reduce the uncertainties related to estimated hydrologic and geophysical properties and predicted hydraulic processes (Hinnell et al., 2010;Mboh et al., 2012).
In the MCMC inversion, it is usually assumed that the noise is additive and Gaussian, leading to the following likelihood function, where N d is the total number of data points, f is a function denoting the forward modeling, and d i and ε i denote the ith data point and associated standard deviation, respectively. The l 1 -norm is used in calculating the likelihood and it is less sensitive to outliers than the l 2 -norm (Tarantola, 1987). It has been observed that most model parameters in this study [θ r , α, n, log(K s ), m, n a , and C sat × 10 −7 ] have limited ranges. For example, the cementation factor m of most geological materials ranges between 1 and 4 (Friedman, 2005). The prior knowledge of the bounds of model parameters can be easily incorporated in calculating the acceptance probability (Equation 25), expressed as where [a, b] is our predefined range for the model parameter z in Z.

SYNTHETIC EXPERIMENT AND INVERSION
In this section, the abovementioned stochastic coupled inversion is used to analyze the information content of transient hydraulic and SP data obtained from a typical outflow experiment. In particular, we will address the following two Note: Q, cumulative outflow; h, pressure head; Δϕ, streaming potential difference; θ r , residual water content; α, fitting parameter interpreted as the inverse of the air-entry pressure; n, the parameter characterizing the shape of the soil water retention curve; K s , saturated hydraulic conductivity; n a , the Archie saturation exponent; C sat , the coupling coefficient at saturation; m, the porosity exponent.

F I G U R E 2
Schematic of a typical outflow experiment. Pore water pressure is measured at two depths using two tensiometers T 1 and T 2 . Self-potential is measured at two depths using a pair of nonpolarizing electrodes E 1 and E 2 . No flux is allowed at the upper boundary, and a pressure head h b is applied to the initially saturated sample to drain the soil. z 1 and z 2 are different elevations questions: (a) what is the uncertainty of the hydraulic and electrical properties estimated from transient hydrogeophysical data; and (b) to what degree will the addition of SP data reduce the uncertainty of the estimated hydraulic properties? Two synthetic soils are considered here (Table 1) for analysis.

Synthetic outflow experiment
The schematic of a typical outflow experiment is shown in Figure 2. The soil sample is initially saturated, and a pressure head h b is applied to the lower boundary to drain the soil sample. The ambient pressure on the upper boundary is equal to the atmospheric pressure but no flux is allowed. During the drainage process, the cumulative outflow Q is monitored. Electrical potential and pressure head are also monitored at two elevations z 1 and z 2 (Figure 2). This type of outflow experiments can be found, for example, in Allègre et al. To test the accuracy of our forward modeling, we conducted an outflow simulation with a cylindrical soil sample (5 cm in radius and 20 cm in height; Figure 2). The sample is initially saturated and the total head is equal to 20 cm with datum at the base of the soil column. The applied h b is −2 m. The cumulative outflow, pore water pressure head, and electrical potential during drainage were calculated using our numerical code. To facilitate the calculation, h b was assumed to decrease linearly from 20 cm to −2 m in a short period (e.g., 100 s). The simulation was stopped when the cumulative outflow Q was not changing significantly. The following soil parameters were used in the simulation: θ s = 0.3 m 3 m −3 , θ r = 0.03 m 3 m −3 , α = 1.4 m −1 , n = 1.6, K s = 5 × 10 −6 m s −1 , n a = 3.5, C sat = −3.5 × 10 −7 V Pa −1 , m = 2, σ s = 0.002 S m −1 , σ w = 0.1 S m −1 , ρ = 1,000 kg m −3 , and g = 10 m s −2 . The electrical potential at the lower boundary is set as zero. The calculated hydraulic responses (h measured between T 1 and T 2 ) and electrical responses (potential difference Δϕ measured between E 1 and E 2 ) are shown in Figure 3. To validate the accuracy of our calculation, the finite element method software COMSOL Multiphysics 5.6 (COMSOL) was also used to conduct the simulation and the results are shown in Figure 3. In general, the results from our numerical modeling are in good agreement with those from the COMSOL simulation. The related RMSDs of Q, h, and Δϕ are 0.88 cm 3 , 0.23 cm, and 0.0028 mV, respectively. This excellent agreement shows the effectiveness and accuracy of our numerical scheme in solving the coupled water flow and streaming potential problem.

Sample 1: Sand
The first soil sample considered here is a sand with a porosity of 0.4. For the soil water retention curve, the van Genuchten model (Equation 2) is used and relevant parameters are α = 8 m −1 , θ s = 0.4, θ r = 0.01, and n = 5; for the unsaturated hydraulic function, Equation 3 is used with log(K s ) = −4 m s −1 . Archie's law incorporating the surface conduction (Equation 8) is used to model the effective electrical conductivity of F I G U R E 3 Transient hydraulic and electrical responses of a soil in an outflow experiment calculated with our numerical code (data points) and COMSOL Multiphysics (solid lines): (a) cumulative outflow Q; (b) pressure head h at two elevations (circles for z 1 and triangles for z 2 ); and (c) streaming potential difference Δϕ between z 1 and z 2 the sand at various saturations with the following parameters: n a = 2.6, m = 1.6, and σ s = 0.0002 S m −1 . We use Equation 15 with C sat = −2.9 × 10 −7 V Pa −1 to model the streaming potential coupling coefficient of the sand with different degrees of saturation. These model parameters are summarized in Table 2. The sample is initially saturated and total head is 20 cm. In the outflow experiment, the pressure head at the lower boundary h b was decreased from the initial value of 20 cm to −1 m within 30 s. The water pressure and electrical potential at z 1 and z 2 are monitored for 10 4 s as well as the cumulative outflow Q. These transient responses are shown in Figure 4. Two scenarios were considered in our stochastic inversion (Table 1). In Scenario 1, only hydraulic data (i.e., Q and h at two elevations) are used in the inversion. The number of F I G U R E 4 Hydraulic and electrical responses of a synthetic sand sample: (a) cumulative outflow Q; (b) pressure head h at two elevations (circles for z 1 and triangles for z 2 ); and (c) streaming potential difference Δϕ between z 1 and z 2 . Dashed and dotted lines represent the predictions calculated using the mean model parameters of Scenarios 1 and 2, respectively data points in each dataset is 30, and thus the total data points used in the inversion is 90. The measurement errors are set as 4 cm 3 and 0.2 cm respectively for Q and h. In the inversion, the parameter θ s is fixed, considering the fact that the porosity of the sand sample is relatively easy to control or measure in the experiment. Initial values of θ r , α, n, and log(K s ) are 0.015, 6 m −1 , 4, and −4.5 m s −1 , respectively, and the covariances are 0.005 2 , 1.5 2 , 1 2 , and 0.2 2 , respectively (also see Table 2). These initial values are selected as typical values of sandy samples (Schaap & Leij, 2000). These initial covariance values are chosen such that the sampling intervals of the proposal distributions (symmetric Gaussian) are approximately half of the predefined ranges of the model parameters (prior knowledge). The MCMC sampling is terminated after 30,000 runs and the chains of the sampled model parameters are shown Note: In the fifth and sixth columns, the estimated mean values of model parameters are in bold, followed by the associated relative differences (in parentheses); the 95% confidence intervals (CIs) are included in square brackets, and the relative lengths of the 95% CIs are in italic. θ r , residual water content; α, fitting parameter interpreted as the inverse of the air-entry pressure; n, the parameter characterizing the shape of the soil water retention curve; log(K s ), logarithm of saturated hydraulic conductivity; n a , the Archie saturation exponent; C sat , the coupling coefficient at saturation; m, the porosity exponent.

F I G U R E 5
The sampled chains of the hydraulic properties of the synthetic sand sample (scenario 1): (a) θ r , (b) α, (c) log K s , and (d) n in Figure 5. The first 10,000 runs are considered as the burnin period, during which the covariances are not updated. The model parameters from the last 20,000 runs are used to estimate the statistical measures of the posterior distributions. The MCMC estimated posterior distributions of the model parameters are shown in Figure 6 as histograms, which are calculated from the last 20,000 runs after the burn-in period. The mean value and 95% confidence interval (CI) of each parameter are estimated from the histograms and listed in Table 2. We also calculated the relative difference (i.e., absolute dif-ference divided by the true value) between the estimated and true mean values as well as the relative length of the 95% CI (i.e., the length of the 95% CI normalized by the true value) for all the model parameters. In Table 2, it is apparent that the recovered model parameters are very close to the true values and the relative differences are generally less than ∼1% except for θ r , which also has a high uncertainty with a relative length of the 95% CI of ∼61.5%. The model parameter log(K s ) has the lowest uncertainty with the length of the 95% CI as ∼0.8% of its mean value. The results in Table 2 Figure 4. It is clear that the simulated Q and h agree very well with "measured" Q and h, and the related RMSDs are only 3.82 cm 3 and 0.19 cm. To evaluate the related prediction uncertainty, a couple of hundred sets of model parameters are drawn randomly from the posterior distributions in Figure 6 and are used to simulate the hydraulic responses. The variation ranges of the simulated Q and h are very narrow (not shown in Figure 4), generally smaller than the size of the data point symbols in Figure 4. Comparisons of experimental and simulated results in Figure 4 confirm that the transient outflow Q and pressure head h data in outflow experiments contain sufficient information on the saturated and unsaturated hydraulic properties of the sample (Toorman et al., 1992).
The second inversion (Scenario 2 in Table 1) is also conducted to estimate both the hydraulic and electrical properties of the sand sample from the "measured" outflow Q, water head h, and streaming potential difference Δϕ. The measurement errors of Q and h are same as those used in Scenario 1, and the error of SP measurements is set as 0.002 mV. In the inversion, the parameters θ s and σ w are considered known; σ s is also assumed known, considering that the effect of surface conduction on SP signal is relatively small . The duration and time intervals of the monitoring data are the same as those used in Scenario 1. Parameters related to the MCMC sampling are also kept unchanged. The prior information on model parameters is summarized in Table 2 as well as the inversion results (mean values and 95% CIs). The inversion results are also shown in Figure 7 as histograms. The estimated mean values are used to simulate the hydraulic and geophysical responses during the drainage and the results are plotted in Figure 4 (dotted lines). Similar to Scenario 1, the agreement between the simulated and "measured" responses is excellent, and the related RMSDs for Q, ∆ϕ, and h are only 4.74 cm 3 , 0.0026 mV, and 0.2 cm, respectively. We also evaluated the variation ranges of the predicted Q, ∆ϕ, and h responses (not shown in Figure 4); compared with Scenario 1, the calculated ranges of variations are quite similar though slightly narrower.
As shown in Table 2, both the recovered hydraulic and electrical parameters in Scenario 2 are very close to the true values. The associated uncertainties (i.e., relative length of the 95% CI) vary between 0.5% for log(K s ) and 30% for θ r . Compared with Scenario 1, the estimated mean values in Scenario 2 are closer to the true values (Table 2); the related uncertainty (relative length of 95% CI) is also smaller, about half of those in Scenario 1. For example, the estimated mean n in Scenario 2 is 5.03 with the 95% CI as [4.97, 5.08], and the estimation from Scenario 1 is 4.97 with the 95% CI as [4.83, 5.11]; the relative length of 95% CI decreases from 5.6 to 2.2%. This indicates that the addition of transient SP data helps the estimation of hydraulic properties. Among the three electrical properties, the cementation factor m has the highest uncertainty, and the relative length of the 95% CI is ∼10.6%; in contrast, the values for C sat and n a are only ∼1 and 1.5%, respectively. This is consistent with the results shown in Younes et al. (2018). The relatively high uncertainty of m is understandable because, in such an outflow experiment (Figure 2), the effective electrical conductivity of the soil is not directly measured. On the contrary, the monitored SP signals are directly influenced by the coupling coefficient C sat . Thus, m is less sensitive to the monitored transient data than C sat (Younes et al., 2018).

Sample 2: Loam
We also consider a loam sample to study if the soil texture affects the MCMC inversion results. Similarly, two scenarios are considered (i.e., inverting with Q and h, and with Q, h, and Δϕ; Table 1). Similarly, θ s , σ w , and σ s are assumed known and the following constant values are used in the inversions: θ s = 0.3 m 3 m −3 , σ w = 0.1 S m −1 , and σ s = 0.002 S m −1 . Similar to the sand sample, the loam sample is initially saturated and the total head is 20 cm. In the forward modeling, a pressure head h b = −2 m is gradually applied to the bottom within  Figure 8). Other parameters are the same as those used in the sand sample. The MCMC inversion results of the loam sample are summarized in Table 3. The transient hydraulic and electrical responses calculated using the MCMC estimated model parameters are shown in Figure 8. Similar conclusions can be made for the loam sample after analyzing the results: (a) using Q and h data are sufficient to accurately estimate the hydraulic properties of the sample; (b) addition of SP data can help reduce the uncertainty of the estimated hydraulic properties of the soil; and (c) the electrical properties of the sample can be estimated from Q, h, and Δϕ data with a fair accuracy, although not as well as hydraulic properties. In general, the lengths of the 95% CIs of the estimated hydraulic properties in Scenario 2 are about one-to two-thirds of those obtained in Scenario 1. Among all the hydraulic properties, θ r has the highest uncertainty, and the relative length of the 95% CI is 96.6% for Scenario 1 and 62.1% for Scenario 2; log(K s ) has the lowest uncertainty, and the relative length of the 95% CI is lower than 2% for both scenarios. For the electrical properties, the estimation of C sat is relatively accurate with the relative length of the 95% CI lower than 4%. In contrast, m and n a have relatively high uncertainties and the relative length of the 95% CIs is higher than 50% (Table 3). Tables 2 and 3 indicates that the estimated model parameters are in general more accurate for the sand than the loam. For example, the relative length of the 95% CI for α is only 2.2% in Scenario 1 and 1.1% in Scenario 2 for the sand, but the values are 7.9 and 5% for the loam. The predicted electrical and hydraulic responses using the recovered model parameters (Figures 4 and 8) show similar trends, particularly for t larger than 100 s. The relatively better performance in sand is probably due to the fact that the transient data of sand used in inversions contain a larger portion of the drainage process (Kool et al., 1985). Note that the last portion of the cumulative outflow Q appears to approach a constant value for the sand (Figure 4), but it still has a tendency to increase for the loam (Figure 8). This implies that, to ensure a good performance on finer soils, the outflow experiment ( Figure 2) should be performed for a longer period than coarse samples such that more information can be included in the measurement.

INFLUENCE OF ELECTRICAL MODELS AND PARAMETERS
It has been found that the unsaturated streaming potential coupling coefficient C model (Equation 15) does not work for many soils (Jougnot et al., 2012Zhang et al., 2017).  In the fifth and sixth columns, the estimated mean values of model parameters are in bold, followed by the associated relative differences (in parentheses); the 95% confidence intervals (CIs) are included in square brackets, and the relative lengths of the 95% CIs are in italic.
θ r , residual water content; α, fitting parameter interpreted as the inverse of the air-entry pressure; n, the parameter characterizing the shape of the soil water retention curve; log(K s ), logarithm of saturated hydraulic conductivity; n a , the Archie saturation exponent; C sat , the coupling coefficient at saturation; m, the porosity exponent.

F I G U R E 8
Hydraulic and electrical responses of the synthetic loam sample: (a) cumulative outflow Q; (b) pressure head h at two elevations (circles for z 1 and triangles for z 2 ); and (c) streaming potential difference Δϕ between z 1 and z 2 . Dashed and dotted lines represent the predictions calculated using the mean model parameters of Scenarios 1 and 2, respectively. The shaded regions represent the variation ranges of the prediction (light color for Scenario 1 and dark color for Scenario 2) In addition, it is also difficult to accurately determine the surface conductivity σ s (Equation 8) of a soil in practice. Thus, it is necessary to analyze how the unsaturated C models and an incorrect σ s affect the coupled inversion.

Influence of an incorrect σ s
We invert the hydraulic and SP measurements of the loam sample once again with σ s = 0, and all the other parameters are kept unchanged. The inverted mean values of the model parameters are θ r = 0.024 m 3 m −3 , α = 1.37 m −1 , n = 1.59, K s = 10 −5.32 m s −1 , n a = 2.33, C sat = −3.51 Note: In the third and fifth columns, the estimated mean values of model parameters are in bold, followed by the associated relative differences (in parentheses); the 95% confidence intervals (CIs) are included in square brackets, and the relative lengths of the 95% CIs are in italic. θ r , residual water content; α, fitting parameter interpreted as the inverse of the air-entry pressure; n, the parameter characterizing the shape of the soil water retention curve; log(K s ), logarithm of saturated hydraulic conductivity; n a , the Archie saturation exponent; C sat , the coupling coefficient at saturation; m, the porosity exponent. ×10 −7 V Pa −1 , and m = 1.63. Compared with the inversion results with an accurate σ s (Scenario 2 in Table 3), most of the recovered parameters are almost unchanged except n a and m, which become noticeably smaller. According to Equation 8, a smaller n a (or m) will increase the calculated electrical conductivity σ of the soil, which will compensate the effect of setting σ s = 0. This implies that using an incorrect σ s in the coupled inversion only affects the other parameters appearing in the electrical conductivity model (Equation 8), and other model parameters are not influenced significantly.

Influence of C models
In this subsection, we analyze if a different coupling coefficient model (other than Equation 15) will affect the conclusion made in our previous section. The unsaturated model proposed in Soldi et al. (2020) is used here, and it can be expressed as (see Equations 23 and 27 in Soldi et al., 2020) = sat e σ sat σ The outflow experiments are simulated for both the sand and loam samples with the new C model and the measurements are then inverted to estimate the petrophysical properties of the samples. In the modeling and inversion, all the model parameters are kept unchanged. The inversion results are shown in Table 4.
In general, the recovered model parameters in Table 4 are very close to those in Tables 2 and 3. Both hydraulic properties and electrical properties of the samples are determined with a very small error. For the sand, the relative differences between estimated and true values range between 0 and 2.8%; for the loam, the related relative differences are between 0 and 10%. Compared with the results of Scenario 1 (Tables 2   and 3), the incorporation of SP data in the inversion with the new C model (Equation 29) also reduces the length of the 95% CIs significantly. For instance, the uncertainty of the parameter α (i.e., relative length of 95% CI) of the loam decreases from 7.9% in Scenario 1 in Table 3 to 2.9% in Table 4. Generally speaking, the inversion results associated with the new C model (Table 4) are consistent with those obtained with Equation 15. Therefore, the conclusions made in this study are still valid if other C models are used in the coupled inversion.

INVERSION OF EXPERIMENTAL DATA
In this section, we apply the stochastic, coupled inversion to experimental data collected from an outflow test reported in Linde et al. (2007). The available datasets include transient Q, h, and ϕ measurements. The soil sample was cylindrical with a radius of 0.035 m and a height of 1.35 m. The sand was initially saturated and the saturated water content θ s (or porosity) is between 0.33 and 0.35 . A pressure head h b = 0.091 m was applied to the bottom to drain the sand column. During the drainage process, both h and ϕ were measured at a number of locations along the sand column. Some hydraulic and electrical properties of the sand have been independently measured  summarized in Table 5): σ sat = 0.012 S m −1 , σ w = 0.051 S m −1 , F = 4.26, C sat = −2.9×10 −7 V Pa −1 , K s = 6.93 × 10 −5 m s −1 .
The datasets used in the inversion include ∼1.8 h of cumulative outflow Q, ∼4 h of pressure head h at two locations (22.5 and 47.5 cm from bottom), and ∼4 h of SP ϕ at two locations (25 and 55 cm from bottom). In total, 30 data points were digitized from Linde et al. (2007) for Q and 40 for h and ϕ; the average time intervals of these data are between 200 and 1,000 s. Similar to synthetic soils, two inversion scenarios were considered (Table 5). In the inversions, the following model parameters were assumed known and thus kept unchanged: θ r = 0.015 m 3 m −3 , σ s = 0, and σ w = 0.051 S m −1 . A total number of 8,000 runs were performed for the MCMC inversions, of which the first 4,000 runs were the burn-in period. The prior information on the model parameters is summarized in Table 5. In Scenario 1, hydraulic parameters were inverted from Q and h and the results are summarized in Table 5. In general, the MCMC estimated mean values are close to the values reported in Linde et al. (2007). For α, n, and log(K s ), the relative differences are less than ∼10%; for θ s , the relative difference exceeds 20%. These discrepancies are probably due to the fact that h and ϕ data used in our inversion only cover ∼4 h of the measurement. In contrast, the inversion in Linde et al. (2007) used transient h and ϕ data with a longer period that includes a free drainage stage (i.e., h b = 9.1 cm was removed) in addition to the initial ∼4 h. Despite these discrepancies, the modeled hydraulic responses using the MCMC estimated model parameters are very close to the measured responses ( Figure 9). The associated RMSDs are only 0.9 cm 3 and 0.55 cm for Q and h, respectively. The good agreement in Figure 9 shows that the stochastic coupled inversion developed in this study could reliably estimate the hydraulic parameters of soils from transient hydraulic data collected in an outflow test.
The inversion results of Scenario 2 are also summarized in Table 5. Comparison of Scenarios 1 and 2 show that the estimated mean values of the hydraulic properties are very similar. The associated uncertainty (e.g., relative length of the 95% CI) is at the same level for parameters θ s , α, n, and log(K s ). These results indicate that the addition of SP data may, although only slightly, reduce the uncertainty of the estimated hydraulic properties of the sand sample.
Regarding the electrical properties, it appears that the estimated mean values of n a , C sat , and m are close to the values reported in Linde et al. (2007), and the related relative differences are 9.4, 17.9, and 4.5%. Compared with hydraulic properties (e.g., results of Scenario 1), the discrepancies of electrical parameters between our estimations and reference values are slightly higher. It is also noted that the estimated 95% CIs (i.e., uncertainty) of the electrical parameters are quite large for n a and m. Similar to synthetic soils, C sat has the lowest uncertainty and the length of the 95% CI is only ∼8.8% of the measured C sat ; this low uncertainty is because SP responses are directly affected by the streaming potential coupling coefficient of the soil (e.g., Equation 10). Moreover, the estimated m value (1.28 in Table 5) is lower than the typical value of unconsolidated sediments (Friedman, 2005). This is because, in the inversion, the surface conductivity was assumed zero (i.e., σ s = 0; Linde et al., 2007). Assuming σ s = 0 will overestimate the contribution of bulk water, thereby underestimating F I G U R E 9 Hydraulic and electrical responses of the sand sample in Linde et al. (2007): (a) cumulative outflow Q; (b) pressure head h at two locations (circles for z 1 = 47.5 cm and triangles for z 2 = 22.5 cm); and (c) self-potential ϕ at two locations (circles for z 1 = 55 cm and triangles z 2 = 25 cm). Dashed and dotted lines represent the predictions calculated using the mean model parameters of Scenarios 1 and 2, respectively the formation factor F and thus the value of m .

CONCLUSIONS
In this study, stochastic, coupled inversions were performed to estimate hydraulic and electrical properties of soils from transient hydrogeophysical data collected in one-step outflow experiments. The time-series data used in the inversion include cumulative outflow Q, pressure head h, and/or SP ϕ. It is found that the stochastic inversion could provide information about the sensitivity of a soil's hydraulic and elec-trical properties to the measured Q, h, and ϕ in the outflow experiment. The results show that both saturated and unsaturated hydraulic properties of soils such as log(K s ), n, and α can be reliably estimated from transient Q and h data and the relative differences between estimation and true values are generally less than 10%. The accuracy of the estimated residual water content θ r is much lower than other hydraulic properties, and the estimation is also associated with a larger uncertainty (the relative length of the 95% CI can reach ∼100%). Comparison of the sand and loam samples indicates that the related uncertainty of the estimated hydraulic parameters can be reduced if the transient hydraulic data used in inversion cover a large portion of the drainage process.
Inversions including SP data were also performed. Compared with inversions with hydraulic data only (i.e., Q and h), the addition of SP data ϕ can reduce the uncertainty of the estimated hydraulic properties, although the mean value is not improved significantly. The related length of the 95% CIs for most hydraulic properties [e.g., log(K s ) and n] could be reduced by about one-third. This indicates that incorporating SP measurement in an outflow experiment is beneficial to the estimation of hydraulic properties. To the best of our knowledge, such a quantitative evaluation of the SP measurement in an outflow experiment was not available in the literature.
Inversion results also show that both saturated and unsaturated electrical properties (C sat , n a , and m) of the soil may be estimated from the hydraulic and SP data collected in outflow experiments. Among these electrical properties, the saturated streaming potential coupling coefficient C sat has the highest accuracy and lowest uncertainty. The accuracy of other parameters such as m and n a is slightly lower compared with C sat (and hydraulic parameters). The associated uncertainty (relative length of the 95% CI) is generally large, about one order of magnitude higher than C sat and hydraulic parameters (e.g., α and n). These different performances are due to the fact that parameters m and n a do not directly affect the generation of streaming potential in an outflow experiment. In addition, it was also found that a good estimation of m and n relies on a correct surface conductivity value in the coupled inversion. The numerical findings in this study can serve as guidance in estimating soil petrophysical properties and their uncertainty from outflow experiments.

A C K N O W L E D G M E N T S
Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for partial support of this research. We thank three reviewers (Damien Jougnot and two anonymous reviewers) and the associate editor (Fred Zhang) for their constructive comments on an earlier version of this paper.