Advancing quantitative understanding of self-potential signatures in the critical zone through long-term monitoring
Kaiyan Hu, Damien Jougnot, Qinghua Huang, Majken C. Looms, Niklas Linde
1 Advancing quantitative understanding of self-potential signatures in the critical zone through long-term monitoring
Kaiyan Hu , Damien Jougnot , Qinghua Huang , Majken C. Looms , Niklas Linde Department of Geophysics, Peking University, 100871 Beijing, China Institute of Earth Sciences, University of Lausanne, 1015 Lausanne Switzerland Sorbonne Université, CNRS, EPHE, UMR 7619 Metis, F-75005, Paris, France Department of Geosciences and Natural Resource Management, University of Copenhagen, Øster Voldgade 10, 1350 Copenhagen, Denmark * Corresponding author. E-mail address: [email protected]
Author contact:
Kaiyan Hu: [email protected]; [email protected] Damien Jougnot: [email protected] Qinghua Huang: [email protected] Majken C. Looms: [email protected] Niklas Linde: [email protected]
This paper has been accepted for publication in Journal of Hydrology:
K. Hu, D. Jougnot, Q. Huang, M. C. Looms, N. Linde (2020), Advancing quantitative understanding of self-potential signatures in the critical zone through long-term monitoring,
Journal of Hydrolog y, 585, 124771, doi: 10.1016/j.jhydrol.2020.124771 ABSTRACT
The self-potential (SP) method is a passive geophysical technique, which may offer insights about water and ionic fluxes in the vadose zone. The main obstacles presently prohibiting its routine use in quantitative vadose zone hydrology are the superposition of signals arising from various source mechanisms, difficult-to-predict electrode polarization effects that depend on electrode design and age, as well as water saturation, pore water chemistry, clay content, and temperature in the immediate vicinity of the electrodes. We present a unique long-term SP monitoring experiment focusing on the first four years of data acquired at different depths in the vadose zone within the HOBE hydrological observatory in Denmark. Using state-of-the-art SP theory combined with flow and transport simulations, we attempt to replicate the observed data and suggest reasons for observed discrepancies. The predictions are overall satisfactory during the first six months of monitoring after which both the patterns and magnitudes of the observed data change drastically. Our main observations are (1) that predicted SP magnitudes are strongly sensitive to how the effective excess charge (or alternatively, the voltage coupling coefficient) scales with water saturation implying that continued research is needed to build more accurate models of electrokinetic phenomena in unsaturated conditions, (2) that significant changes in electrode polarization occur in the shallowest electrodes at time scales of a year, most likely due to desaturation by capillarity of the fluids filling the electrodes, suggesting that electrode effects cannot be ignored and that explicit electrode modeling should be considered in future monitoring studies, and (3) that multi-rate mass transfer and reactive transport modeling, with specific emphasis on fluid-mineral interactions, are needed to better predict salinity and pore water conductivity. We hope to stimulate other researchers to test new SP modeling approaches and interpretation strategies against these data by making the SP and complimentary (temperature, dielectric constant, potential/actual evapotranspiration, precipitation) data time-series available. Keywords:
Hydrogeophysics; Long-term monitoring; Self-potential; Vadose zone; Electrode effects; Solute transport. Introduction
Self-potential (SP) signals are mainly generated by naturally-occurring electrokinetic, electrochemical and bioelectrical source mechanisms (Revil and Jardani, 2013). In the SP method, passive measurements with non-polarizing electrodes are used to map or monitor voltages on the ground surface or in boreholes. Due to a rich set of natural source mechanisms, the SP method has been used to study a wide range of hydrological and critical zone processes. Most studies have based their interpretations on the electrokinetic effect, which originates when water flowing in porous media drags excess electric charges in the electrical double layer to produce the so-called streaming potential. SP signals attributed to streaming potentials have been used to characterize hydraulic parameters (e.g., Zlotnicki and Nishida, 2003; Jardani et al., 2009; Jougnot et al., 2012; Soueid Ahmed et al., 2014), landslides (Lapenna et al., 2003) and seepage from dams (e.g., Ikard et al., 2014; Rittgers et al., 2014; Soueid Ahmed et al., 2019). SP signals of electrochemical origin arise in the presence of ionic concentration gradients in the pore water or when redox reactions are facilitated by electronic conductors. Self-potential signals of electrochemical origin have been used to constrain saline tracer and contaminant transport (e.g., Maineult et al., 2005; Martínez-Pagán et al., 2010; Jougnot et al., 2015), sea water intrusion (e.g., Graham et al., 2018; MacAllister et al., 2018), to locate the corrosion of buried metallic bodies (e.g., Castermant et al., 2008; Rittgers et al., 2013) and to invert for redox maps associated with contaminant plumes (Linde and Revil, 2007). In the SP literature, it is often tacitly ignored that recorded voltages do in situ water fluxes in the vadose zone? Or more simply, can SP data be used to determine if water fluxes at a given time and depth are directed upward or downward? We contrast the measured SP data with state-of-the-art SP numerical modeling that accounts for electrokinetic and electrochemical contributions, as well as electrode effects. In an effort to engage the SP community in research towards an improved mechanistic and quantitative understanding of SP signals in the critical zone, we make the full SP and complimentary (dielectric constant, temperature, potential and actual evapotranspiration, precipitation) data record available. Methodology
Section 2.1 introduces basic self-potential theory followed by more detailed sections describing the electrokinetic contribution (section 2.1.1), the electro-diffusive contribution (section 2.1.2), and electrode-effects (section 2.1.3) leading to the total SP response (section 2.1.4). Section 2.2 introduces the hydrodynamic characterization and modeling that is needed input for self-potential modeling in the vadose zone. Water flow and solute transport (section 2.2.1) is first introduced before explaining how in situ field data are used to infer vadose zone properties by inversion (section 2.2.2).
Below we describe the theory used to perform numerical simulations of SP signals in the vadose zone. In this treatment, we exclude the redox contribution since we do not expect electronic conductors at the HOBE site. In the absence of external current sources, we have : ∇ ∙ 𝐉 𝐭𝐨𝐭 = ∇ ∙ (𝐉 𝐜 + 𝐉 𝐞𝐤 + 𝐉 𝐞𝐝 + 𝐉 𝐞𝐞 ) = 0 , (1) where the total electric current density 𝐉 𝐭𝐨𝐭 is the sum of the conductive current density 𝐉 𝐜 (A/m ) , the electrokinetic current density 𝐉 𝐞𝐤 (A/m ) and the electro-diffusive current density 𝐉 𝐞𝐝 (A/m ) and electrode-related contributions 𝐉 𝐞𝐞 (A/m ) . The conductive current density 𝐉 𝐜 (A/m ) obeys Ohm’s law: 𝜎(𝜃 w )𝐄 = 𝐉 𝐜 , (2) where 𝐄 (V/m) is the electric field, which is the negative gradient of electrical potential 𝜑 (V), 𝜎(𝜃 w ) (S/m) is the effective electrical conductivity of the soil. For coarser soils without fines, we have (Archie, 1942) : 𝜎(𝜃 w ) = 𝐹 −1 𝜎 w 𝑆 𝑤𝑛 , (3) where 𝑆 𝑤 = 𝜃 w 𝜙 is the water saturation, 𝜃 w is the water content, 𝐹 = 𝜙 −𝑚 is the formation factor, 𝜙 is the porosity, 𝑚 is the cementation exponent and 𝑛 is the saturation exponent. For soils containing fines, the relationship between bulk electrical conductivity and pore water electrical conductivity is rather complex and non-linear (e.g., Revil et al., 2019). Assuming that the pore water is a NaCl-solution, the electrical conductivity of pore water 𝜎 w (S/m) , the temperature of water 𝑇 (℃) and solute concentration 𝐶 w (mol/L) are related by (Sen and Goode, 1992): 𝜎 w = (5.6 + 0.27𝑇 − 1.5 × 10 𝑇 )𝐶 w − w 𝐶 w . (4) The electrokinetic current density depends on the water flux 𝐮 and the effective excess charge density 𝑄 veff ̂ (𝜃 w , 𝐶 w ) (C/m ) of porous media: 𝐉 ek = 𝑄 veff ̂ (𝜃 w , 𝐶 w )𝐮 , (5) where 𝑄 veff ̂ (𝜃 w , 𝐶 w ) varies, for a given soil, with the water content 𝜃 𝑤 and solute concentration 𝐶 w . It might also exhibit a hysteretic behavior that we ignore in the present treatment (Revil et al., 2008). The main challenge in deriving the effective excess charge is which it is a flux-averaged property (Linde, 2009). In this work, we will mainly rely on the water retention (WR) formalism introduced by Jougnot et al. (2012) in which the soil is conceptualized as a bundle of capillaries that can be either filled with water or air (Jackson, 2008). In this approach, 𝑄 WReff ̂(𝜃 w , 𝐶 w ) is derived from the water retention function by integrating the effective excess charge density 𝑄 veff,𝑟 ̂ (𝑟, 𝐶 w ) of each saturated capillary with radius 𝑟 over all capillaries in the soil: 𝑄 WReff ̂(𝜃 w , 𝐶 w ) = ∫ 𝑄 veff,r ̂ (𝑟,𝐶 w ) rmaxrmin 𝑣 𝑟 (𝑟)𝑓 WR (𝑟)𝑑𝑟∫ 𝑣 𝑟 (𝑟) rmaxrmin 𝑓 WR (𝑟)𝑑𝑟 , (6) where 𝑣 𝑟 (𝑟) is the pore water flux in a capillary and 𝑓 WR (𝑟) is the capillary size distribution derived from the relationship between pore pressure and saturation under the van Genuchten water retention model (van Genuchten, 1980). An alternative for calculating the capillary size distribution is the relative permeability (RP) method 𝑓 RP (𝑟) , which derives the saturation from the relative permeability function (Jougnot et al., 2012). For the WR-approach, the relative excess charge density 𝑄 WReff,rel ̂ (𝜃 w , 𝐶 w ) curves defined by the ratio of effective excess charge density with saturated excess charge density 𝑄 WReff ̂ (𝜃 w ,𝐶 w )𝑄 WReff,sat ̂ (𝐶 w ) depend on the shape of the water retention curves ( 𝛼 VG and 𝑛 VG ) of different materials. This approach has been shown to predict the evolution of the effective excess charge density with respect to the water content (e.g., Jougnot and Linde, 2013; Voytek et al., 2019). Here, 𝛼 VG ( m −1 ) refers to the inverse of the air entry pressure and 𝑛 VG (-) is related to the pore-size distribution. Instead of using Eq. (6) and in perfect analogy with the treatment of hydraulic conductivity in unsaturated media, we scale this relative excess charge density with laboratory or field data at saturated conditions. This leads to 𝑄 eff ̂ (𝜃 w , 𝐶 w ) = 𝑄 WReff,rel ̂ (𝜃 w , 𝐶 w ) 𝑄 eff,sat ̂ (𝐶 w ) , (7) where 𝑄 eff ̂ (𝜃 w , 𝐶 w ) are the values used in subsequent modeling. Guarracino and Jougnot (2018) and Soldi et al. (2019) proposed analytical solutions to obtain these parameters, but with overall similar results. A common form of the empirical relationship between the voltage coupling coefficient 𝐶 eksat (mV/m) at saturation and electrical conductivity of the pore water 𝜎 w (S/m) is (Revil et al., 2003; Linde et al., 2007b; Jougnot et al., 2015): log |𝐶 eksat | = 𝑎 + 𝑏 log 𝜎 w + (𝑐 log 𝜎 w ) , (8) where Jougnot et al. (2015) found that 𝑎 = −0.6 , 𝑏 = −1.319 , 𝑐 = −0.1227 are appropriate for modeling SP profiles at the HOBE hydrological site. From this, 𝑄 eff,sat ̂ (𝐶 w ) needed in Eq. (7) is obtained as (Revil and Leroy, 2004) : 𝑄 eff,sat ̂ = − 𝐶 eksat (V/m)𝜎𝜌 w 𝑔𝐾 sat , (9) where electrical conductivity σ (S/m) is computed by Eq. (4), 𝜌 w (kg/ m ) is the density of groundwater, 𝑔 (m/ s ) is the acceleration of gravity and 𝐾 sat (m/s) is the saturated hydraulic conductivity. Differences in ion mobility 𝛽 𝑗 lead to electro-diffusive currents in the presence of electrochemical potential gradients. Considering Na + and Cl − as the only ions (𝐶 w =𝐶 Na + = 𝐶 Cl − ) in the pore water, the electro-diffusive current density is written as (Pride, 1994; Revil and Jougnot, 2008; Jougnot et al., 2015): 𝐉 𝐞𝐝 = − k B 𝑇 K e (𝑡 NaH − 𝑡
ClH )𝜎(𝜃 w )∇ln (𝐶 w ) , (11) where 𝑇 K is the temperature in Kelvin (K), e = 1.602 × 10 −19 (C) is the elementary charge and k B = 1.3806 × 10 −23 (J/K) is the Boltzmann constant . Furthermore, 𝑡 NaH and 𝑡 ClH are the microscopic Hittorf transport numbers related to the ion mobility 𝛽 Na and 𝛽 Cl : 𝑡 NaH = 𝛽 Na 𝛽 Na +𝛽 Cl , 𝑡 ClH = 1 − 𝑡
NaH . (12) The microscopic Hittorf transport number 𝑡 NaH depends on concentration at high salinity, while it is constant at low salinities (Revil, 1999; Gulamali et al., 2011; Leinov and Jackson, 2014): 𝑡 NaH = { 0.3962, 𝐶 w < 0.09 mol/L0.3655 − 9.2 × 10 −3 ln(𝐶 w ) , 𝐶 w ≥ 0.09 mol/L . (13) When considering microporous media (e.g., clay), macroscopic Hittorf transport numbers (Revil and Jougnot, 2008) should be used to properly account for the anion exclusion effect (Leinov and Jackson, 2014). At the HOBE site, according to the investigation of grain-size distribution (Uglebjerg, 2013), the soil mostly consists of sand (> 90%), and the clay barely occurs in the shallow 50 cm (< 1%) and in the range of 100-250 cm depth (0.5% to 2%). We can, therefore, rely on the microscopic Hittorf transport numbers assuming a negligible cation exchange capacity (CEC) (see Jougnot et al. 2015 for comparative discussion of alternative formulations). Non-polarizable electrodes are electrodes for which the electrode polarization does not depend on the current flowing through the electrode. At the HOBE site, we rely on Petiau electrodes (Pb/PbCl , NaCl) that are described in detail by Petiau (2000) and manufactured by SDEC (PMS9000). These are electrodes of the second kind, which implies that an auxiliary salt (NaCl or KCl) is used in addition to the principal salt (PbCl ). Petiau electrodes have been designed to ensure long-term stability (the electrode drift is on the order of 0.2 mV/month), minimal noise and low temperature dependence. The choice of Pb/PbCl as a metal-salt couple is motivated by an extensive comparison with other alternatives by Petiau and Dupis (1980). The stability of Petiau electrodes is partly due to an appropriate combination of pH, and high principal salt (PbCl ) and auxiliary salt (here NaCl, with KCl being a popular alternative) concentrations. Furthermore, the outlet in the form of a porous plug is located comparatively far (approximately 0.15 m) from the lead wire, and the design includes an internal thin channel located behind the porous plug enabling electrolytic contact while decreasing ionic diffusion. To further decrease ionic diffusion, the electrolyte is mixed with clay mud. With this design, the desaturation time, defined as the time for the diluted front originating at the outlet to reach the lead wire, is expected to be on the order of 14 years, assuming saturated conditions. Experiments mentioned by Petiau (2000) suggests that this time might be further increased by 20-40% when electrodes are installed under partially saturated conditions. This suggests that we can assume that the ionic concentrations of the principal and auxiliary salts in the vicinity of the lead wire are constant over many years. Nevertheless, this does not imply that the potential difference between the lead wire and the pore water outside the porous plug is constant over time. A standard consideration is that the electrode potential depends on temperature, which implies that temperature differences between electrodes lead to unwanted voltages that need to be removed. For small temperature ranges, the self-potential due to temperature differences can be expressed by a linear relationship (Petiau, 2000): 𝑆𝑃 temp = 𝑘 T (𝑇 obs − 𝑇 ref ), (14) where the temperature coefficient 𝑘 T depends on the electrode type. For Pb/PbCl Petiau electrodes (Petiau, 2000) with a saturated NaCl electrolyte, 𝑘 T is , 𝑇 obs is the temperature in the position of the measurement electrode and 𝑇 ref is the temperature at the reference electrode. To predict the temperature at any depth and time, we use the measured temperature data to fit a spatial-temporal distribution of temperature with a sine function: 𝑇(𝑡, 𝑧) = 𝐴(𝑧) sin{𝜔[𝑡 − 𝜑 (𝑧)]} + 𝑇 c , (15) where 𝜔 = (rad/d) is the angular frequency, 𝑇 c is the annual average temperature that is assumed constant. The amplitude factor 𝐴(𝑧) and the phase delay 𝜑 (𝑧) is considered to vary linearly with depth. To further improve the corrections particularly at shallow depths, we add the interpolated residual values 𝑟𝑒𝑠(𝑡, 𝑧) as: 𝑇̃ (𝑡, 𝑧) = 𝑇(𝑡, 𝑧) + 𝑟𝑒𝑠(𝑡, 𝑧) . (16) A second non-standard consideration is that the membrane potential between the interior of the clay mud at concentration C ele and the pore water outside of the porous pot ( C obs and C ref for the concentrations in the vicinity of the observation and reference electrodes, respectively) depends on the pore water chemistry. For practical purposes, we are interested in the differential membrane potential 𝑆𝑃 mb , that is, the difference in membrane potential (Revil, 1999) between a measuring and reference electrode: 𝑆𝑃 mb = − k B 𝑇 K e (2𝑇 Na∗ − 1)ln 𝐶 ele 𝐶 obs − [− k B 𝑇 K e (2𝑇 Na∗ − 1)ln 𝐶 ele 𝐶 ref ] = k B 𝑇 K e (1 − 2𝑇 Na∗ )ln 𝐶 ref 𝐶 obs , (17) where 𝑇 Na∗ is the macroscopic Hittorf transport number in the clay mud filling the Petiau electrodes. Determining this value is difficult as the detailed design of the mud mixture in the PMS9000 electrodes is unknown. Using Petiau (2000) as a guide, we make the following assumptions concerning the clay mud in the electrodes: cation exchange capacity CEC = 2.89 × 10 C/kg (0.03 in meq/g), density 𝜌 e = 2600 kg/m , porosity 𝜙 e = 0.5 and cementation exponent 𝑚 e = 1.5. With these values, the total excess of charge 𝑄 v can be calculated by (Waxman and Smits, 1968): 𝑄 v = 𝜌 e 1−𝜙 e 𝜙 e CEC . (18) This volumetric excess charge 𝑄 v of the clay mud is needed to determine the surface electrical conductivity 𝜎 s and is different in nature than the flux-averaged effective excess charge density described in Eq. (5) (see discussions in Jougnot et al., 2019; Jougnot et al., 2020). The ratio of 𝜎 s to the electrical conductivity of electrolyte 𝜎 w is (Sen and Kan, 1987; Revil et al., 1998): 𝛾 e = 𝜎 s 𝜎 w ≈ s 𝑄 v w ( 𝜙 e e ) = s ρ e CEC3𝜎 w , (19) where 𝛽 s is the clay surface ionic mobility taken as −9 (m s −1 V −1 ) for Na + (Revil, 1999). The macroscopic Hittorf number 𝑇 Na∗ grows as a function of this ratio 𝛾 e following (Revil, 1999): 𝑇 Na∗ = [ 1 + Na∗ 𝛾 e 𝜙 e−𝑚e + (𝑡 Na∗ −𝛾 e )(1− 𝛾e𝑡Na∗ +√(1− 𝛾e𝑡Na∗ ) + ) ] −1 , (20) where the microscopic Hittorf number 𝑡 Na∗ of Na + in the clay mixture of the electrode can be estimated by Eq. (13) assuming that 𝐶 e is 6.8 mol/L (Petiau, 2000). Based on the above assumptions, 𝑇 Na∗ of the electrode clay mixture is 0.3492, which is only slightly larger than 𝑡 Na∗ due to the high concentration.
The electrical potential response of the streaming current 𝐉 ek and the electro-diffusive source 𝐉 ed are computed and indicated as 𝑆𝑃 ek and 𝑆𝑃 ed , respectively. The total SP signal is: 𝑆𝑃 = 𝑆𝑃 ek + 𝑆𝑃 ed + 𝑆𝑃 temp + 𝑆𝑃 mb . (21) We design 𝑆𝑃 sim = 𝑆𝑃 ek + 𝑆𝑃 ed + 𝑆𝑃 mb as the simulated SP response that is to be compared with the temperature-corrected raw data using Eqs. (14) and (15) that we refer to as 𝑆𝑃 obs . A modified version of the finite difference MaFlot code (Künze and Lunati, 2012) is used to simulate the SP response by solving the Poisson’s equation obtained by combining Eqs. (1) and (2). Figure 1a displays a sketch to illustrate the calculation of total SP. To apply the framework of section 2.1, we need inputs in terms of distributed time-series of water content, water fluxes and water concentrations. To achieve this, we use a finite element solution of the Richards equation and the advection-dispersion equation. Richards equation is used to simulate water flow under variably-saturated conditions (Richards, 1931): 𝜕𝜃 w 𝜕𝑡 + 𝑆 = ∇(𝐾∇𝐻 p ) , (22) where 𝑆 ( s −1 ) is the source term, 𝐻 p (m) is the water pressure head, and 𝐾 (m/s) is the hydraulic conductivity expressed as: 𝐾(𝜃 w ) = 𝐾 rel (𝜃 w )𝐾 sat , (23) where 𝐾 rel (-) is the relative hydraulic conductivity described by the van Genuchten (1980) model with parameters 𝛼 VG and 𝑛 VG . For given boundary and initial conditions, 𝐻 p and 𝜃 w can be calculated by the governing equation Eq. (22). The resulting water flux 𝐮 is expressed as: 𝐮 = − 𝐾(𝜃 w )∇𝐻 p . (24) For solute transport, we rely on the advection-dispersion equation (ADE): 𝜕(𝜃 w 𝐶 𝑗 )𝜕𝑡 + ∇ ∙ [𝐮𝐶 𝑗 − 𝜃 w (𝐷 𝑗 + 𝛼v̅)∇𝐶 𝑗 ] = 0 , (25) where v̅ = 𝐮𝜃 w is the average flow velocity, 𝐷 𝑗 ( m /s ) is the ionic diffusion coefficient and 𝛼 (m) is the dispersivity (Bear, 2012). Using a conservative ADE, the fresh infiltrating rainwater will in a recharge-dominated setting eventually fill the saturated portion of the pore space. This is unrealistic as demonstrated by the fact that the electrical salinity of the undisturbed pore water is around 0.02 S/m (Jougnot et al., 2015). Coupling the flow- and transport simulations with a predictive soil chemical model including interactions between the soil minerals and the water phase is extremely challenging and outside the scope of this study. In an effort, to account for the leading effect of how fresh rainwater is turned into more charged pore water, we use a very simple iterative approach (Fig. 1a). We output the profile of simulated concentrations 𝐶 cur daily ( 𝑑𝑡 = 1 day ) and compare with an assumed equilibrium concentration 𝐶 eq : if 𝐶 cur < 𝐶 eq then 𝐶 in = 𝛾𝑑𝑡(𝐶 eq − 𝐶 cur ) + 𝐶 cur , (26) else 𝐶 in = 𝐶 cur , where 𝛾 determines how quickly equilibrium is reached. Considering two extreme conditions: if 𝛾 = 0 d −1 , no adaptation is made and 𝐶 in = 𝐶 cur ; if 𝛾 = 1 d −1 , we have 𝐶 in = 𝐶 eq . The resulting 𝐶 in is used as the initial concentration when simulating the following day. We stress that this adaptation in the modeling framework is made to qualitatively assess the impact of rain water properties on SP phenomena in the near surface and that it is not a consistent physics-based model. Based on Linde et al. (2006), we express the effective dielectric constant 𝜅 𝑡sim as 𝜅 𝑡sim = (𝐹 − 1)𝜅 s + (𝜃 w 𝜙 ) 𝑛 𝜅 w + [1 − (𝜃 w 𝜙 ) 𝑛 ] 𝜅 a 𝐹 , (27) where 𝜅 s = 4.3 and 𝜅 a = 1 are the dielectric constants of the solid grains and air, respectively. 𝜅 w is the dielectric constant of pore water, which depends on temperature 𝑇 (℃) following (Weast et al., 1988; Wraith and Or, 1999): 𝜅 w = 78.54[1 − 4.579 × 10 −3 (𝑇 − 25) + 1.19 × 10 −5 (𝑇 −25) − 2.8 × 10 −8 (𝑇 − 25) ] . (28) The observed temperature data at the field site is used in Eq. (28) to estimate the dielectric constant of water 𝜅 w before using Eq. (27) to obtain the simulated effective dielectric constant 𝜅 𝑡sim from the simulated water contents. We contrast the 𝜅 𝑡sim with the observed effective dielectric constant 𝜅 𝑡obs to evalute the likelihood of a given vadose zone model by using an ℓ -norm of residual errors. This leads to the following Laplacian log-likelihood function: ℒ( 𝜃 wr , 𝜙, 𝛼 VG , 𝑛 VG , 𝐾 sat |𝜅 obs ) = − ∑ ∑ (ln (2𝜎 𝑡𝑖 )) 𝑀𝑖=1𝑁𝑡=1 − (29) ∑ ∑ ( |𝜅 𝑡obs,𝑖 −𝜅 𝑡sim,𝑖 |𝜎 𝑡𝑖 ) 𝑀𝑖=1𝑁𝑡=1 , where 𝜎 𝑡𝑖 is a data error, M is the total number of available sensors and 𝑁 is the total number of times considered. The log-prior density is also taken into account when calculating the log-posterior density of each model realization. 𝑝(𝜃 wr , 𝜙, 𝛼 VG , 𝑛 VG , 𝐾 sat |𝜅 obs ) = ℒ(𝜃 𝑤𝑟 , 𝜙, 𝛼 𝑉𝐺 , 𝑛 𝑉𝐺 , 𝐾 sat |𝜅 obs ) +∑ ∑ ln[ℕ(𝜇 𝑖𝑘 , 𝜎 𝑖𝑘 , 𝑋 𝑖𝑘 )] , (30) where for q layers ℕ(𝜇 𝑖𝑘 , σ 𝑖𝑘 , 𝑋 𝑖𝑘 ) represents the normally-distributed prior distribution with 𝜇 𝑖𝑘 and σ 𝑖𝑘 the expected value and standard deviation of parameter 𝑘 at the 𝑖th layer, 𝑋 𝑖𝑘 represents the sampled parameter value of 𝜃 wr , 𝜙, 𝛼 VG , 𝑛 VG and 𝐾 sat respectively. The inversion process follows the diagram in Fig. 1b. As indicated by Fig. 1b, we apply Markov chain Monte Carlo (MCMC) to infer the maximum a posteriori estimate of the model parameters. For this, we rely on the DREAM (ZS) algorithm (Laloy and Vrugt, 2012; Vrugt, 2016). We refer to Jougnot et al. (2015) for more details. Fig. 1. Flow diagrams for simulating SP signals (modified from Jougnot et al., 2015) displaying (a) the SP simulation scheme and (b) the inversion scheme for hydraulic parameters.
3. The HOBE observatory and data
The HOBE hydrological observatory is located within an agricultural field (Fig. 2) in the Skjern river catchment (Fig. 2a), near Voulund, western Denmark (Jensen and Illangasekare, 2011; Jensen and Refsgaard, 2018). It was established to study hydrological responses generated by natural and artificial forcings at different scales and to establish water balance closure in the context of assessing and managing water resources within the catchment by integration of data from multiple measuring techniques (e.g., geophysics, hydrology, meteorology, remote sensing, geochemistry). The patterns of dynamic hydrological landsurface and subsurface processes are analyzed based on multi-scale monitoring in space and time. Previous findings relevant to this work include (1) a description of the geological architecture and soil components based on analysis of drilling samples (Uglebjerg, 2013); (2) the estimation of potential and actual evapotranspiration (Ringgaard et al., 2011) and the establishment of the water balance at the field-scale (Vásquez et al., 2015); (3) the detection of soil moisture and other pools of hydrogen using the cosmic-ray neutron method (Andreasen et al., 2016; Andreasen et al., 2017); (4) monitoring the pathway of saline movement based on cross-borehole electrical resistivity tomography (ERT) and ground-penetrating radar (GPR) (Haarder et al., 2015) and from SP data (Jougnot et al., 2015). Based on the analyzed samples from five drill cores (Uglebjerg, 2013), the soil is described by seven horizontal layers. In each layer, more than 90% of the soil is made up of very fine to very coarse sand, with silt and clay content in layers 1 and 2 being 5-10%. In the present work, we assume that surface conductivity is negligible given the soil description and that we did not have access to induced polarisation data (e.g., Revil et al., 2019) at the site. The water table is fluctuating between 5.5 m and 6.5 m depth. In September 14, 2011, 40 liters of saline water with an electrical conductivity of 22.1-24 S/m was injected during two hours uniformly over a 12 m × 12 m experimental plot. The tracer experiment was described and analyzed by Haarder et al. (2015) and the SP signature of the saline injection was studied by Jougnot et al. (2015). Fig. 2. (a) Map of Denmark with (b) top view with locations of sensors and (c) depth information about the sensors (modified from Jougnot et al., 2015). Electrode SP-0.5B was installed several meters away from the area depicted in (b).
Precipitation data have been collected hourly at the site since April 2009 using multiple rain gauges. A time series was developed, combining data from five rain gauges, in order to obtain what was named the
Best Assembled Precipitation Dataset . In this dataset, the majority of the data (~70%) came from a ground-based rain gauge located within 25 m from the SP installations (see below). Hourly potential evapotranspiration data have been calculated with the Penman-Monteith equation (Monteith, 1965) using wind speed, temperature, humidity, solar radiation and crop height data from November 2008 until August 2013. Based on the analysis of rain water, Uglebjerg (2013) found that its average electrical conductivity is approximately 50 μS/m . In order to predict the SP response, we first use a similar approach (see section 2.2.2) as Jougnot et al. (2015) to calibrate the hydraulic parameters, with the main difference that we rely on a longer time-record (December 1, 2011 to July 22, 2013) of dielectric constant and temperature data. We kept the petrophysical parameters 𝑚 and 𝑛 values in Eq. (27) to those obtained by Jougnot et al. (2015), namely 1.38 and 1.57, respectively. For each tested combination of soil parameter values, Hydrus-1D (Šimůnek et al., 2013) was used to calculate the water content 𝜃 w at the sensor locations (0.25 m, 1 m, 1.5 m, 2 m, 2.25 m, 2.5 m, 2.75 m, and 3 m). We initialize our simulations on August 1, 2010 assuming a pressure head profile of -50 cm. The time-varying flux boundary condition at the surface is determined by precipitation and potential evapotranspiration data, except for the period of saline injection in September 2011 during which the vertical inflow was imposed at 3.6 cm/d. At the bottom of the numerical domain with thickness of 6.5 m, a fixed pressure head of 50 cm is used corresponding to an assumed constant water table at 6 m. The profile is divided into 7 soil layers that are discretized using 300 nodes. From the flow simulations, we output the water content 𝜃 w daily at midnight in the period from December 1, 2011 to July 22, 2013. Table 1:
Uncorrelated prior Gaussian distribution of model parameters (modified from Jougnot et al., 2015) with 𝜇 and 𝜎 indicating the mean value and the standard deviation, respectively. 𝐾 sat is saturated hydraulic conductivity and 𝜃 𝑤𝑟 is residual water content. Layer
Depth ( 𝐦 ) 𝜽 𝐰𝐫 𝝓 𝐥𝐨𝐠 𝜶 𝐕𝐆 (𝐜𝐦 −𝟏 ) 𝒏 𝐕𝐆 𝐥𝐨𝐠 𝑲 𝐬𝐚𝐭 (𝐜𝐦/𝐝) 𝜇 𝜎 𝜇 𝜎 𝜇 𝜎 𝜇 𝜎 𝜇 𝜎
1 0-0.2 0.05 0.01 0.39 0.01 -0.87 0.5 1.36 0.1 3.31 1.47 2 0.2-0.45 0.04 0.01 0.38 0.02 -1.15 0.32 2.3 0.36 2.74 0.49 3 0.45-0.8 0.03 0.01 0.38 0.02 -1.29 0.1 2.3 0.36 2.44 0.46 4 0.8-1.2 0.07 0.03 0.4 0.02 -1.02 0.17 1.71 0.31 2.13 0.12 5 1.2-1.75 0.06 0.02 0.37 0.04 -0.96 0.14 2.49 0.64 2.46 0.43 6 1.75-2.55 0.04 0.01 0.39 0.02 -1.21 0.36 2.89 1.04 2.46 0.43 7 2.55-7.3 0.04 0.01 0.39 0.02 -1.29 0.22 2.89 1.04 2.42 0.43 Gaussian prior distributions of material properties were used (Table 1) and the initial models were obtained as random draws from the prior. The Metropolis acceptance rule (Metropolis et al., 1953) is used to determine, at each MCMC step, if the proposal is accepted. We use three MCMC chains each having a length of 2000. The number of layers q is 7 in Eq. (30). The standard deviation of the dielectric constant data, 𝜎 𝑡𝑖 , appearing in Eq. (29) data is assumed constant and equal to 1. This value is not so crucial as we are interested in the maximum a posteriori estimate; it was determined by trial-and-error and describes how well we are able to fit our observed data. The log-likelihood value increases rapidly during the first 200 MCMC steps before increasing more gradually. Using Eq. (30), we calculate the corresponding log-posterior densities and the maximum value is picked. The corresponding model parameter values are then used in the subsequent SP simulations (Table 2). Figure 3 is a scatter plot comparing the simulated and measured dielectric constant data for this best-fitting model. The average absolute mean deviation is 0.7168. The simulated dielectric constant underestimates the variability in the observed data, particularly during the winter season when the water content can be very high due to decreased evapotranspiration and snow melt events (Haarder et al., 2015). Table 2:
The parameter values of the maximum a posteriori model obtained by MCMC inversion.
Layer
Depth ( 𝐦 ) 𝜽 𝒘𝐫 𝝓 𝜶 𝐕𝐆 (𝐜𝐦 −𝟏 ) 𝒏 𝐕𝐆 𝑲 𝐬𝐚𝐭 (𝐜𝐦/𝐝)
1 0-0.2 0.066
2 0.2-0.45 0.071
3 0.45-0.8 0.023
4 0.8-1.2 0.020
5 1.2-1.75 0.064
6 1.75-2.55 0.069
7 2.55-7.3 0.053 Fig. 3. Scatter plot of simulated dielectric constant against measured dielectric constant in the (a) winter season (October to March) and the (b) summer season (April to September) of 2012.
With assuming only vertical flow, a total of 15 non-polarizable Petiau electrodes (Pb/Pb Cl , NaCl ) were installed at depths of 0.25 m, 0.5 m, 0.75 m, 1 m, 1.45 m, 1.9 m, 3.1 m and 7.3 m in individual holes (Fig. 2c). Theses electrodes holes were drilled by augering. Each electrode was inserted vertically in a thin hole and it was attached to a wooden stick to ensure that the electrode had a good contact with the bottom of the hole and to better estimate the depth of installation. The one exception to this procedure was one of the top electrodes that was installed at 0.5 m in an individual dug hole, which was later filled with bentonite. As shown in Fig. 2c, there are two electrodes installed without bentonite surroundings at 0.25 m, 0.5 m, 0.75 m, 1 m and 1.45 m. The reference electrode during the measurements was located well below the water table at 7.3 m depth. To remove the strong influence of water table fluctuations from the acquired SP data, they were processed to have the deepest electrode in the vadose zone (at 3.1 m depth) serve as the reference electrode. Due to the superposition principle, this is achieved by subtracting the measured SP data at this location from all the other measured SP time-series. Initially, from July 19, 2011 until September 9, 2011, SP data were recorded every 20 minutes. Subsequently, SP data were recorded every 5 minutes. Herein, we present data until July 28, 2015. Repeated crosshole ERT measurements were performed in the area from September 2011 to August 2012 (Haarder et al., 2015), leading to large recorded voltages that are unrelated to the processes of interest. To filter out these effects, we first median filter the SP data at each recorded time 𝑡 𝑖 within a period of 𝑡 𝑖 − 75 minutes to 𝑡 𝑖 + 75 minutes. Second, in the period of ERT measurements we identify periods ∆𝑡 𝑘 disturbed by the ERT measurements by identifying times with anomalously high temporal SP variations. For these time-periods, we replace the raw data with the median filtered data. Finally, we median filter the resulting time-series using a 1-hour window corresponding to the time-interval used to represent the SP data in this paper. See Mariethoz et al. (2015) for an alternative data processing approach that has been demonstrated for a sub-set of the presented data. Results
We first consider the streaming potential generated by electrokinetic effects. In doing so, we focus on the SP data acquired prior to the saline tracer injection. In this period immediately following the installation, the electrodes are expected to be stable and the electro-diffusive effect should be minimal. In addition to simulations using the WR and RP methods mentioned in section 2.1.1, we also consider modeling based on the thick double-layer assumption in which the effective excess charge density is scaled with water saturation 𝜙 𝜃 w (Linde et al., 2007a): 𝑄 Lindeeff ̂ (𝜃 w , 𝐶 w ) = 𝜙 𝜃 w 𝑄 eff,sat ̂ (𝐶 w ) , (31) where the saturated excess charge density 𝑄 eff,sat ̂ (𝐶 w ) is estimated by Eqs. (8) and (9). We present SP profiles (Figs. 4b-d) at three times corresponding to before, during and after a rainfall event (Fig. 4a). As expected for a recharge-dominated regime with associated downward flow, the recorded SP data (Figs. 4b-d) become increasingly negative towards the surface. Furthermore, the electrodes located at the same depth display similar SP values indicating that electrode-related effects are uniform. The simulated SP signals have similar shapes (Figs. 4b-d), but it is only the WR method that provides simulations of the same magnitude as the observed data. The RP method overestimates the amplitudes roughly by a factor of two and the thick double-layer assumption of Eq. (31) leads to simulated SP signals that are roughly 50% too small. The temporal differences in the observed SP data (Figs. 4e-g) display small decreases in SP magnitudes, which is in contrast to the simulations that show larger increases. In the following, we only present predictions based on the WR method. In the simulations described above, we used calibrated values of the inverse time constant 𝛾 = 0.025 d −1 (Eq. 26) obtained by maximizing the correlation between simulated and observed SP signals. As shown in Figs. 4h-j, the vertical SP profiles are only slightly affected by the delay factor for 𝛾 > 0.025 d −1 . Fig. 4. (a) Precipitation data with three chosen times (t1-t3) highlighted by black arrows. Comparison between observed and simulated SP profiles (b) before a heavy rain at t1, (c) during the rain at t2 and (d) after the rain at t3. The solid lines are simulated SP data using different excess charge densities 𝑄 eff ̂ obtained, respectively, from the WR method (red), the RP method (blue) and the thick double-layer model (green) of Linde et al. (2007). (e-g) Differences of observed SP (circles) and simulated SP (dashed lines) with respect to t1. (h-j) Simulated SP signals (solid lines) using the WR method with different γ -values in equation (26). After choosing the WR method for calculating the effective charge density (Eq. 7) and calibrating the delay factor at 𝛾 = 0.025 d −1 (Eq. 26) using the observed SP data prior to saline injection, we now consider the SP simulations for longer time-series including the saline tracer injection period. Following the ERT-based estimations of the saline tracer test by Haarder et al. (2015), the dispersivity α (Eq. 25) is initially set to 0.25 m. We compare the SP simulations with the observed SP signals from July 2011 to September 2012 at 0.25 m and 0.75 m (Fig. 5). The simulated SP data given as a superposition of electrokinetic and electro-diffusive contributions without electrode-related effects are presented in Figs. 5b and 5e. The simulated SP signals describe the main shape of the observed SP data rather well, particularly the strong response caused by the saline tracer injection. To assess electrode-related effects, we first consider the estimated macroscopic Hittorf number 𝑇 Na∗ of the electrode clay-mixture calculated by Eq. (20). The simulated SP signals including the electrode-effect (Figs. 5c and 5f) display a much slower decrease in SP amplitudes after the tracer injection than what is observed, and the gradual relaxation observed in the data towards more negative values is not reproduced. We then sought the macroscopic Hittorf number 𝑇 Na∗ that maximises the correlation between the observed and simulated SP data. In Figs. 5d and 5g, we display the prediction for the best-fitting calibrated value of 0.43. This higher value is rather consistent with the fact that the clay allows a larger mobility of the cation due to an important volume of EDL in its pore space (see discussion in Revil and Jougnot, 2008). Here, the simulated responses are found to be intermediate between the two other cases. For the shallow electrode, we see that the decrease of SP magnitudes (Fig. 5d) following the saline tracer injection is sharper than for the other cases (Figs. 5c-f) and in better agreement with the observed data. Fig.5. (a) Precipitation data during the first year of SP monitoring with the red arrow indicating the time of saline tracer injection. Solid black and dashed lines indicate the observed SP data at two measuring electrodes located at (b-d) 0.25 m and (e-g) 0.75 m depth. (b, e) The sum of the simulated streaming potential and electro-diffusive potential is shown in blue. The sum of the simulated streaming potential, electro-diffusive potential and electrode membrane potential with (c, f) the calculated macroscopic Hittorf number of the electrode clay mixture in purple and (d, g) with the best fitting macroscopic Hittorf number (0.43) in red.
Using the optimized 𝑇 Na∗ = 0.43 for the electrode clay mixture, we now investigate the sensitivity of our simulated SP responses to the dispersivity α . The impact of α at 0.25 m (Fig. 6a) and 0.75 m (Fig. 6b) mainly concerns the speed at which the SP magnitudes decrease following the tracer injection, with faster decreases observed for larger α , and the rate of relaxation towards more negative values being the strongest for small α . The most satisfactory simulations are offered by α = 0.1 m and this value is retained in the following. As an aside, we note that the observed SP curves at 0.25 m (Fig. 6a) present a notable decrease in February 2012, which is attributed to the partially frozen topsoil (and electrode). Indeed, the soil temperature at this depth is below 0 ℃ in this time period. In Fig. 7, we present the simulated spatio-temporal evolution of the solute concentration, which suggests that all of the injected saline should have reached the groundwater in April 2013. The trajectory of the simulated maximum concentration agrees well with the ERT-based estimates (Haarder et al., 2015), with both showing an acceleration of the plume movement in the winter. We note further that our estimates are in-between the ERT-based estimates of Haarder et al. (2015) and those determined by analyzing core data. Fig. 6. Comparison between the observed (black solid and dashed lines) and simulated SP signals using the best fitting macroscopic Hittorf number (0.43) for different values of dispersivity at (a) 0.25 m and (b) 0.75 m depth. The chosen values of dispersivity are 0.05 m (green), 0.1 m (purple), 0.15 m (blue) and 0.25 m (red). The red arrow indicates the time of saline tracer injection and the blue arrows indicate the frozen period. Fig. 7. The spatio-temporal variations of the simulated solute concentration of pore water assuming a dispersivity of 0.1 m. The dashed line corresponds to the depth of maximum concentration, the black triangles indicate the center of solute mass determined by drilled core data and the red diamonds are the center of solute mass inferred by Haarder et al. (2015) using ERT data.
Using the calibrated parameters, we simulate the SP responses as a superposition of electrokinetic, electro-diffusive, and electrode-related effects over four years and consider five vertical SP profiles in the summer seasons of 2011, 2012, 2013, 2014, and 2015 (Fig. 8). The corresponding dates were chosen to correspond to periods with heavy rainfall implying that the simulated SP data have similar shapes and magnitudes at all of the considered dates. However, it is only in 2011 (Fig. 8a) shortly after the SP installations that the observed SP data also show a trend of increasingly negative values towards the surface as expected for an infiltration-dominated system. At later years, the observed data start to deviate from this expected behavior. The observed data are negative in 2012 (Fig. 8b), but the amplitudes are higher and there are rather significant oscillations with depth. This pattern is similar in 2013 (Fig. 8c), but with smaller magnitudes. The observed SP data turn positive in 2014 (Fig. 8d) and the magnitudes of the most positive values in 2015 approach 200 mV (Fig. 8e), thus, the different scale used on x-axis of Fig. 8e. Since no major changes in the state of the vadose zone is expected between these years as indicated by the simulated data, this suggests that the observed differences between the years are largely due to time-varying electrode polarizations manifesting themselves by electrode-potentials that develop differently with depth (different water content magnitudes and dynamics). Fig. 8. Comparison between observed (black dots) and simulated (red lines) SP profiles on (a) 07 September 2011, (b) 22 August 2012, (c) 10 August 2013, (d) 30 August 2014 and (e) 28 July 2015 (e), respectively. These dates in during the late summer period of each year were chosen to correspond to days just during important rain events (25.32 mm, 33.77 mm, 18.66 mm, 22.48 mm and 37.74 mm, respectively).
Based on the simulation of solute transport presented in Fig. 7, the injected saline water should disappear from the vadose zone after April 2013. Figure 8 highlighted strong deviations from the simulated data starting a year from the electrode installation that we are unable to explain without invoking changes in the electrode potential (i.e., electrode drifts). The common approach when facing unexplained electrode offsets, is to refer to changes in the SP data with respect to a reference datum. This approach can only work satisfactorily if the “electrode drift” is slow and unrelated to the local water content and salinity. With the aim of highlighting the SP response to heavy rainfall events after more than three years of monitoring, we consider the SP data at 0.25 m and 1.9 m during three days in the summer of 2014 and the winter of 2015, respectively. We show the curves in Figs. 9 and 10 with different scales, as there are large differences in amplitudes between 0.25 m and 1.9 m. In Figs. 9 and 10, the SP data are presented as relative changes compared with the SP data recorded at the corresponding depth at 00:00:00 on August 30, 2014 and at 12:00:00 on January 7, 2015, respectively. This type of re-referencing of SP data is very common (e.g., Linde et al., 2007a). Furthermore, the reference electrode was set at 2.5 m to focus on temporal changes in the near surface. The periods were chosen to correspond to a long dry period followed by heavy rainfall (Fig. 9a and Fig. 10a) leading to strong changes in water content (Figs. 9b, d and Figs. 10b, d). In the summer season of 2014, one of the SP electrodes at 0.25 m show responses corresponding to the timing of changes in the water content. At 1.9 m, there is an SP response (Fig. 9e) with the same sign as in the winter (Fig. 10e), but the initiation of the SP signal occurs with a time delay with respect to the water content signal (Fig. 9d). Considering the winter period, the SP data at 0.25 m (Fig. 10c) do not show a clear unambiguous relationship to the rainfall event and the two SP electrodes respond rather differently. However, the SP data at 1.9 m (Fig. 10e) show as expected increasing negative amplitudes starting at a time that is only slightly earlier than the response in the water content at 2.0 m (Fig. 10d). These differences might suggest significant lateral variations in vertical flow. Fig. 9. SP responses observed during a period with significant rain events separated by dry periods in the summer of 2014. (a) Precipitation data; (b) measured water content data at 0.25 m; (c) observed SP data at 0.25 m; (d) measured water content data at 2 m and (e) observed SP data at 1.9 m. Here, the SP reference is located at 2.5 m depth and the amplitudes are set to zero at the beginning of the plotted period. Fig. 10. SP responses observed during a period with significant rain events separated by dry periods in the winter of 2015. (a) Precipitation data; (b) measured water content data at 0.25 m; (c) observed SP data at 0.25 m; (d) measured water content data at 2 m and (e) observed SP data at 1.9 m. Here, the SP reference is located at 2.5 m depth and the amplitudes are set to zero at the beginning of the plotted period. Discussion
The electrode surrounded by bentonite was installed several meters away from the saline injection area and is likely unaffected by the saline injection. The bentonite with its low permeability is likely to offer a stable environment in which water content and salinity variations are low. Indeed, the use of bentonite is often recommended for SP monitoring (Petiau, 2000). We consider here the SP data at 0.5 m, from the installation time until the saline tracer injection. In this time period, it is expected that the main signals are of electrokinetic origin. The overall shape of the SP electrode in bentonite (Fig. 11b) is similar to the numerical simulations, which ignores the presence of the bentonite (Fig. 11d). The two time-series show a stepped shape with declines occurring around July 22, August 7 and September 6 corresponding to the main rain events in this time-period. On top of this, the electrode in the bentonite shows a shift in magnitude and a linear drift with increasing negative values over time that is not present in the simulated data. One of the electrodes installed at 0.5 m (Fig. 11c) shows temporal variations corresponding to the same events, but the magnitudes are smaller than in the bentonite and the polarity of the response is not always the same (negative on July 22, negative on August 7 and positive on September 6). Fig. 11. (a) Precipitation data in the considered time period. Comparison between the observed SP signals in (b) bentonite with (c) observed SP signals without bentonite installation and (d) simulated SP ignoring the presence of the bentonite at 0.5 m. SP data were not recorded from 26 August 2011 to 30 August 2011.
A notable feature is clearly that the SP signals in the bentonite (Fig. 11b) are much smaller (-65 mV on average) than the other SP electrodes at the same depth (Fig. 11c) and the simulated SP data (Fig. 11d). This is likely related to the fact that bentonite is a microporous media with a much larger fraction of the pore space being affected by the electric double layers (EDL). One would thus expect anion exclusion effects (Jougnot et al., 2009; Leinov and Jackson, 2014) implying that Cl − is excluded from the pore space, leading to a different macroscopic Hittorf number 𝑇 Naben in the bentonite compared to the surrounding. Considering the simple case for which the soil salinity is uniform, the difference of SP between the bentonite 𝑆𝑃 benexc and the surrounding soil 𝑆𝑃 soilele is: ∆𝜑 benexc = 𝑆𝑃 benexc − 𝑆𝑃 soilele = − k B 𝑇 K e (2𝑇 Na∗ − 1)ln 𝐶 ele 𝐶 ben − k B 𝑇 K e (2𝑇 Naben − 1)ln 𝐶 ben 𝐶 soil + k B 𝑇 K e (2𝑇 Na∗ − 1)ln 𝐶 ele 𝐶 soil = k B 𝑇 K e (1 − 2𝑇 Na∗ )ln 𝐶 soil 𝐶 ben + k B 𝑇 K e (1 − 2𝑇 Naben )ln 𝐶 ben 𝐶 soil (32) where 𝐶 ele , 𝐶 ben and 𝐶 soil are the ion concentrations inside the electrode, in the bentonite and in the surrounding soil. Considering the limiting case 𝑇 Naben = 1 , ∆𝜑 benexc could be -66.32 mV when 𝐶 ben 𝐶 soil = 10 and 𝑇 Na∗ = 0.43 at 20 ℃ (293.15 K) in good agreement with the observations in Figs. 11b and 11c. This example suggests that the electrode installation method has a large effect on the measured SP signals (absolute amplitudes, drifts and sensitivity to perturbations) and that monitoring of ion concentrations or suitable proxies (e.g., electrical conductivity) in the vicinity of the electrodes are needed. One risk in the context of multi-year SP monitoring is that other ions than Cl − that are naturally present in the soil enter into contact with the lead wire through diffusion (Petiau, 2000). This may lead to unpredictable perturbations of the electrode polarization that arise well before the desaturation time of Cl − (Petiau, 2000). Judging from our data (Fig. 8), an even larger risk is related to dehydration of the electrolyte by capillarity occurring when the SP electrodes are installed under unsaturated conditions with the dehydration process varying over time and space as the capillary pressure varies. In contrast to the short-term laboratory findings by Petiau (2000) discussed in section 2.1.3, our long-term field-based results suggest that the lifespan of Petiau electrodes is greatly reduced in partially saturated media. We attribute this to capillary suction removing the electrolyte from the clay mud inside the electrodes. Petiau (2000) reports that the measured electrode potential increases when the electrodes dry out as a consequence of the limited contact between the remaining electrolyte and the lead wire inside the electrodes. This is in strong agreement with the observed shift towards positive values over time (Fig. 8). Our field observations suggest that Petiau electrodes installed in the vadose zone have considerably shorter life-times than previously suggested by Petiau (2000). It could be argued that the impact of time-varying saturation and concentration state on the electrodes could be included in the modeling and that their impact on the measured SP data could be predicted. However, observed differences between electrodes located close to each other at the same depth are essentially unpredictable when considering 1-D flow and transport simulations. To characterize this aspect of electrode drift related to local variations in the electrode installation, geological heterogeneity and perhaps the electrodes themselves, we plot the median of the absolute deviation of the electrode pairs installed at common depths (Fig. 12a). The electrodes are stable with low absolute deviations before May 2012, after which the discrepancies increase and reach a peak of almost 40 mV around May 2013. It is also seen that the median of the absolute deviations is overall lower in the summers than in the winters. The largest absolute deviations are found for the shallow electrodes. This behavior is also seen in Fig. 5, in which the deviations among the electrode pairs increase notably after May 2012. The median value of the absolute deviation between simulated water content and measured water content inverted by measured dielectric constant data using Eq. (27) shows similar seasonal trends (Fig. 12b). Figure 12b offers, thus, some possible insights regarding the possible causes of the larger absolute deviations in the observed SP data (Fig. 12a) during the winter period. Indeed, partial freezing of the soil and snow accumulation on the surface lead to heterogeneous melt water infiltration (making the underlying 1-D assumption increasingly questionable) and various measurement issues (Vásquez et al., 2015). Fig. 12. (a) Temporal variation of the absolute median deviation of observed SP data at depths (0.25 m, 0.5 m, 0.75 m, 1 m, 1.45 m) with two measurement electrodes. The corresponding absolute median deviations of each electrode pair for the whole time period are 15.08 mV, 12.92 mV, 9.03 mV, 9.04 mV, 8.96 mV, respectively. (b) Temporal variation of the absolute median deviation of observed water content data with simulated water content data at depths (0.25 m, 1 m, 1.5 m, 2 m, 2.25 m, 2.5 m, 2.75 m, 3 m). The corresponding absolute median deviations at each depth over time are 0.03, 0.02, 0.03, 0.01, 0.02, 0.02, 0.02, 0.02, respectively. Concluding remarks
The analysis of four years of SP monitoring data at the HOBE hydrological observatory in Denmark has highlighted that modeling of streaming potentials in the vadose zone is still a largely unresolved challenge. The main complication is the strong variability of the effective excess charge ( 𝑄 eff ̂ (𝜃 w , 𝐶 w ) ) and the difficulties in predicting this function accurately. Traditionally, it is expected that streaming potential magnitudes increase when water fluxes increase (e.g., Doussan et al., 2002). However, this does not need to be the case and increased flow might actually lead to decreases in the observed SP amplitudes (see Figs. 4f and 4g). As the water content increases, 𝑄 eff ̂ (𝜃 w , 𝐶 w ) decreases particularly at the onset of fluid flow in macro pores (these pores carry a lot of flow without a significant drag of excess charge). So, this decrease in 𝑄 eff ̂ (𝜃 w , 𝐶 w ) might offset the increase in current density offered by the increasing water flux u (Eq. 5), while the increasing water content increases electrical conductivity (Eq. 3), thereby, decreasing the observed self-potential (Eq. 2). Calibration of such non-linear functions as 𝑄 eff ̂ (𝜃 w , 𝐶 w ) to achieve high predictive capability is very challenging. Dual-domain formulations of porosity need to be considered to account for the influence of macro pores on flow, transport and electrical properties (see also discussion Romero-Ruiz et al. 2019). This work has also highlighted the need to accurately model the time-evolving water salinity, which implies that water-soil interactions need to be accounted for, and that water salinity needs to be monitored. Furthermore, it is probably also needed to account for the fact that flux-averaged soil water concentrations are different from in-situ concentrations (e.g., Evaristo et al., 2015). This implies that predictive SP modeling in the vadose zone does not only need to consider advanced SP theory, but also state-of-the-art flow and transport modeling accounting for chemical soil interactions in a dual- or multi-continuum framework. Our work clearly highlights that SP signals in the vadose zone cannot be interpreted as a sensor of water movement by simply differencing SP differences between neighboring depth levels. The reality is much more complicated than this even in absence of electrode effects. Our study confirms suggestions by Perrier et al. (1997) and Doussan et al. (2002) that the design and installation of SP-electrodes for long-term monitoring in the vadose zone is a largely open topic. Important electrode-related challenges in predictive SP modeling concern ( i ) the membrane polarization between the interior of the electrode and its immediate surroundings and ( ii ) dehydration of the electrodes by capillary forces impacting electrode potentials. Except for one electrode, we did not follow the common approach of installing the electrodes in a bentonite solution (Perrier et al., 1997; Petiau, 2000; Doussan et al., 2002). Clearly, a bentonite solution improves the contact resistance with the ground and makes the data less sensitive to small-scale heterogeneity (Petiau, 2000). However, it also leads to difficult-to-predict electro-diffusive effects as the bentonite equilibrates with the surrounding media and the electrode solution leading in our case to large electrode offsets with respect to other electrodes and enhanced drifts (Fig. 11b). For instance, the soil mud in which Doussan et al. (2002) inserted their electrodes displayed a threefold decrease in the pore water electrical conductivity during the monitoring period and, thus, a similar decrease in salinity. This implies significant electro-diffusive effects of 3-D nature in the soil and does not ensure constant water salinity at the contact with the electrode. For future studies, we strongly recommend that pore water conductivity and chemistry is monitored, which would help to elucidate some of the behaviors that we report. Acknowledgments
We thank the Danish hydrological observatory HOBE for providing access to the site, technical help, and full access to the data. The data used in this paper are available here: http://dx.doi.org/10.17632/6r8898657w.1. To calculate the SP signals, we used a version of Maflot (maflot.com) that was kindly provided by I. Lunati and R. Künze and modified by M. Rosas-Carbajal. We are grateful for the financial support offered by the China Scholarship Council and the helpful reviews by Emily Voytek and André Revil.
References
Andreasen, M., Jensen, K. H., Desilets, D., Zreda, M., Bogena, H., Loom, M. C., 2017. Cosmic-ray neutron transport at a forest field site: the sensitivity to various environmental conditions with focus on biomass and canopy interception. Hydrol. Earth Syst. Sci. 21, 1875-1894. https://doi.org/10.5194/hess-21-1875-2017. Andreasen, M., Jensen, K. H., Zreda, M., Desilets, D., Bogena, H., Loom, M. C., 2016. Modeling cosmic ray neutron field measurements. Water Resour. Res. 52, 6451-6471. https://doi.org/10.1002/2015WR018236. Archie, G., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Min. Metall. Eng. 146, 54-61. https://doi.org/10.2118/942054-G. Bear, J., 2012. Hydraulics of groundwater. Courier Corporation. Castermant, J., Mendonça, C. A., Revil, A., Trolard, F., Bourrié, G., Linde, N., 2008. Redox potential distribution inferred from self-potential measurements associated with the corrosion of a burden metallic body. Geophys. Prospect. 56(2), 269-282. https://doi.org/10.1111/j.1365-2478.2007.00675.x. DesRoches, A. J., Butler, K., 2016. Monitoring and modelling of pumping-induced self-potentials for transmissivity estimation within a heterogeneous confined aquifer. Geophys. J. Int. 207, 1722-1738. https://doi.org/10.1093/gji/ggw354. Doussan, C., Jouniaux, L., Thony, J., 2002. Variations of self-potential and unsaturated water flow with time in sandy loam and clay loam soils. J. Hydrol. 267, 173-185. https://doi.org/10.1016/S0022-1694(02)00148-8. Evaristo, J., Jasechko, S., McDonnell, J. J., 2015. Global separation of plant transpiration from groundwater and streamflow. Nature. 525, 91-94. https://doi.org/10.1038/nature14983. Graham, M. T., MacAllister, D., Vinogradov, J., Jackson, M. D., Butler, A. P., 2018. Self-potential as a predictor of seawater intrusion in coastal groundwater boreholes. Water Resour. Res. 54. http://doi.org/10.1029/2018WR022972. Guarracino, L., Jougnot, D., 2018. A physically based analytical model to describe effective excess charge for streaming potential generation in water saturated porous media. J. Geophys. Res. Solid Earth 123(1), 52-65. https://doi.org/10.1002/2017JB014873. Gulamali, M. Y., Leinov, E., Jackson, M. D., 2011. Self-potential anomalies induced by water injection into hydrocarbon reservoirs. Geophysics 76(4), F283-F292. https://doi.org/10.1190/1.3596010. Haarder, E. B., Jensen, K.H., Binley, A., Nielsen, J., Uglebjerg, T. B., Looms, M. C., 2015. Estimation of recharge from long-term monitoring of saline tracer transport using electrical resistivity tomography. Vadose Zone J. 14(7). https://doi.org/10.2136/vzj2014.08.0110. Ikard, S. J., Revil, A., Schmutz, M., Karaoulis, M., Jardani, A., Mooney, M. A., 2014. Characterization of focused seepage through an earthfill dam using geoelectrical methods. Ground Water 52(6), 952-965. https://doi.org/10.1111/gwat.12151. Jackson, M. D., 2008. Characterization of multiphase electrokinetic coupling using a bundle of capillary tubes model. J. Geophys. Res. Solid Earth 113(B04201). https://doi.org/10.1029/2007JB005490. Jardani, A., Revil, A., Barrash, W., Crespy, A., Rizzo, E., Straface, S., Cardiff,M., Malama, B., Miller, C., Johnson, T., 2009. Reconstruction of the water table from self-potential data: a bayesian approach. Ground Water 47(2), 312-227. https://doi.org/10.1111/j.1745-6584.2008.00513.x. Jensen, K. H., Illangasekare, T. H., 2011. HOBE: A Hydrological Observatory. Vadose Zone J. 10(1), 1-7. https://doi.org/10.2136/vzj2011.0006. Jensen, K. H., Refsgaard, J. C., 2018. HOBE: The Danish Hydrological Observatory. Vadose Zone J. 17(1), 1-24. https://doi.org/10.2136/vzj2018.03.0059. Jougnot, D., Linde, N., 2013. Self-potentials in partially saturated media: the importance of explicit modeling of electrode effects. Vadose Zone J. 12(2): 1-49. https://doi.org/10.2136/vzj2018.03.0059. Jougnot, D., Linde, N., Haarder, E. B., Looms, M. C., 2015. Monitoring of saline tracer movement with vertically distributed self-potential measurements at the HOBE agricultural test site, Voulund, Denmark. J. Hydrol. 521, 314-327. https://doi.org/10.1016/j.jhydrol.2014.11.041. Jougnot, D., Linde, N., Revil, A., Doussan, C., 2012. Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils. Vadose Zone J. 11(1), 272-286. https://doi.org/10.2136/vzj2011.0086. Jougnot, D., Mendieta, A., Leroy, P., Maineult, A., 2019. Exploring the effect of the pore size distribution on the streaming potential generation insaturated porous media, insight from pore network simulations. J. Geophys. Res. Solid Earth 124(6), 5315-5335. https://doi.org/10.1029/2018JB017240. Jougnot, D., Revil, A., Leroy, P., 2009. Diffusion of ionic tracers in the Callovo-Oxfordian clay-rock using the Donnan equilibrium model and the formation factor. Geochimica et Cosmochimica Acta, 73(10), 2712-2726. https://doi.org/10.1016/j.gca.2009.01.035. Jougnot, D., Roubinet, D., Guarracino, L., Maineult, A., 2020. Modeling streaming potential in porous and fractured media, description and benefits of the effective excess charge density approach, accepted for publication in the book “Advances in Modeling and Interpretation in Near Surface Geophysics" Eds.: A. Biswas and S.P. Sharma, Springer Geophysics Series, ISBN 978-3-030-28909-6. Künze, R., Lunati, I., 2012. An adaptive multiscale method for density-driven instabilities. J. Comput. Phys. 231, 5557-5570. https://doi.org/10.1016/j.jcp.2012.02.025. Laloy, E., Vrugt, J. A., 2012. High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing. Water Resour. Res. 48(1), W01526. https://doi.org/10.1029/2011WR010608. Lapenna, V. et al., 2003. High-resolution geoelectrical tomographies in the study of Giarrossa landslide (southern Italy). Bull. Eng. Geol. Env. 62(3), 259-268. https://doi.org/10.1007/s10064-002-0184-z. Leinov, E., Jackson, M. D., 2014. Experimental measurements of the SP response to concentration and temperature gradients in sandstones with application to subsurface geophysical monitoring. J. Geophys. Res. Solid Earth 119(9), 6855-6876. https://doi.org/10.1002/2014JB011249. Linde, N., 2009. Comment on ‘‘Characterization of multiphase electrokinetic coupling using a bundle of capillary tubes model’’ by Mathew D. Jackson. J. Geophys. Res. Solid Earth 114, B06209. https://doi.org/10.1029/2008JB005845. Linde, N., Binley, A., Tryggvason, A., Pedersen, L. B., Revil, A., 2006. Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data. Water Resour. Res., 42(11), W12404. https://doi.org/10.1029/2006WR005131. Linde, N., Jougnot, D., Revil, A., Matthai, S. K., Arora, T., Renard, D., Doussan, C., 2007a. Streaming current generation in two-phase flow conditions. Geophys. Res. Lett. 34(3), L03306. https://doi.org/10.1029/2006GL028878. Linde, N., Revil, A., 2007. Inverting self‐potential data for redox potentials of contaminant plumes. Geophys. Res. Lett., 34(14), L14302. https://doi.org/10.1029/2007GL030084. Linde, N. Revil, A., Bolève, A., Dagès, C., Castermant, J., Suski, B., Voltz, M., 2007b. Estimation of the water table throughout a catchment using self-potential and piezometric data in a Bayesian framework. J. Hydrol. 334, 89-99. https://doi.org/10.1016/j.jhydrol.2006.09.027. MacAllister, D., Jackson, M. D., Butler, A. P., Vinogradov, J., 2018. Remote detection of saline intrusion in a coastal aquifer using borehole measurements of self‐potential. Water Resour. Res. 54(3), 1669-1687. https://doi.org/10.1002/2017WR021034 Maineult, A., Bemabé, Y., Ackerer, P., 2005. Detection of advected concentration and pH fronts from self-potential measurements. J. Geophys. Res. Solid Earth, 110, B11205. https://doi.org/10.1029/2005JB003824. Mariethoz, G., Linde, N., Jougnot, D., Rezaee, H., 2015. Feature-preserving interpolation and filtering of environmental time series. Environ. Modell. Softw. 72: 71-76. https://doi.org/10.1016/j.envsoft.2015.07.001. Martínez-Pagán, P., Jardani, A., Revil, A., Haas, A., 2010. Self-potential monitoring of a salt plume. Geophysics 75(4), WA17-WA25. https://doi.org/10.1190/1.3475533. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E., 1953. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys., 21, 1087-1092. https://doi.org/10.1063/1.1699114. Monteith, J. L., 1965. Evaporation and environment. Symp. Soc. Exp. Biol. 19, 205-224. Park, S. K., Larsen, J. C., Lee, T-C., 2007. Electrical resistivity changes not observed with the 28 September 2004 M6.0 Parkfield earthquake on the San Andreas fault, California. J. Geophys. Res. 112, B12305. https://doi.org/10.1029/2006JB004459. Perrier, F. Petiau, G., Clerc, G., Bogorodsky, V., Erkul, E., Jouniaux, L., Lesmes, D., Mancnae, J., Meunier, J., Morgan, D., Nascimento, D., Oettinger, G., Schwarz, G., Toh, H., Valiant, M., Vozoff, K., Yazici-Cakin, O., 1997. A one-year systematic study of electrodes for long period measurements of the electric field in geophysical environments. J. Geomag. Geoelectr. 49, 1677-1696. https://doi.org/10.5636/jgg.49.1677. Petiau, G., 2000. Second generation of lead-lead chloride electrodes for geophysical applications. Pure Appl. Geophys. 157(3), 357-382. https://doi.org/10.1007/s000240050004. Petiau, G., Dupis, A., 1980. Noise, temperature coeefficient, and long time stability of electrodes for telluric observations. Geophys. Prospect. 28, 792-804. https://doi.org/10.1111/j.1365-2478.1980.tb01261.x. Pride, S., 1994. Governing equation for the coupled electromagnetics and acoustics of porous media. Phys. Rev. B 50(21), 678-696. https://doi.org/10.1103/PhysRevB.50.15678. Revil, A., Cathles III, L. M., Losh, S., Nunn, J. A., 1998. Electrical conductivity in shaly sands with geophysical applications. J. Geophys. Res. Solid Earth 103(B10), 23925-23936. https://doi.org/10.1029/98JB02125. Revil, A., 1999. Ionic diffusivity, electrical conductivity, membrane and thermoelectric potentials in colloids and granular porous media: A unified model. J. Colloid Interface Sci. 212(2), 503-522. https://doi.org/10.1006/jcis.1998.6077. Revil, A., Gevaudan, C., Lu, N., Maineult, A., 2008. Hysteresis of the self-potential response associated with harmonic pumping tests. Geophys. Res. Lett. 35, L16402. https://doi.org/10.1029/2008GL035025. Revil, A., Jardani, A., 2013. The self-potential method. Cambridge University Press, New York. Revil, A., Jougnot, D., 2008. Diffusion of ions in unsaturated porous materials. J. Colloid Interface Sci, 319: 226-235. https://doi.org/10.1016/j.jcis.2007.10.041. Revil, A., Leroy, P., 2004. Constitutive equations for ionic transport in porous shales. J. Geophys. Res. Solid Earth 109, B03208. https://doi.org/1029/2003JB002755. Revil, A., Naudet, V., Nourzaret, J., Pessel, M., 2003. Principles of electrography applied to self-potential electrokinetic sources and hydrogeological applications. Water Resour. Res., 39(5), 1114. https://doi.org/10.1029/2001WR000916. Revil A., Qi, Y., Ghorbani, A., Coperey, A., Soueid Ahmed, A., Finizola, A., Ricci, T., 2019. Induced polarization of volcanic rocks. 3. Imaging clay cap properties in geothermal fields. Geophys. J. Int., 218(2), 1398-1427. https://doi.org/10.1093/gji/ggz207. Richards, L. A., 1931. Capillary conduction of liquids through porous media. Physics 1, 318-333. https://doi.org/10.1063/1.1745010. Ringgaard, R. Herbst, M., Friborg, T., Schelde, K., Thomsen, A. G., Soegaard, H., 2011. Energy fluxes above three disparate surfaces in a temperate mesoscale coastal catchment. Vadose Zone J. 10, 54-66. https://doi.org/10.2136/vzj2009.0181. Rittgers, J. B., Revil, A., Karaoulis, M., Mooney, M. A., Atekwana, E. A., 2013. Self-potential signals generated by the corrosion of buried metallic objects with application to contaminant plumes. Geophysics, 78(5), 65-82. https://doi.org/10.1190/GEO2013-0033.1. Rittgers, J. B., Revil, A., Planes, T., Mooney, M. A., Koelewijn, A. R., 2014. 4-D imaging of seepage in earthen embankments with time-lapse inversion of self-potential data constrained by acoustic emissions localization. Geophys. J. Int. 200(2), 758-772. https://doi.org/10.1093/gji/ggu432. Romero-Ruiz, A., Linde, N., Keller, T., Or, D., 2019. A review of geophysical methods for soil structure characterization. Rev. Geophys. 56, 672-697. https://doi.org/10.1029/2018RG000611. Sen, P. N., Goode, P. A., 1992. Influence of temperature on electrical conductivity on shaly sands. Geophysics 57, 89-96. https://doi.org/10.1190/1.1443191. Sen, P. N., Kan, R., 1987. Electrolytic conduction in porous media with charges. Phys. Rev. Lett., 58(8): 778-780. https://doi.org/10.1103/PhysRevLett.58.778. Šimůnek, J., Šejna, M., Saito, H., Sakai, M., van Genuchten, M. T., 2013. The Hydrus-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media. Version 4.17, HYDRUS Software Series 3. Department of Environmental Sciences, University of California Riverside, Riverside, California, USA. Soldi, M., Jougnot, D., Guarracino, L., 2019. An analytical effective excess charge density model to predict the streaming potential generated by unsaturated flow. Geophys. J. Int. 216(1), 380-394. https://doi.org/10.1093/gji/ggy391. Soueid Ahmed, A., Jardani, A., Revil, A., Dupont, J. P., 2014. Hydraulic conductivity field characterization from the joint inversion of hydraulic heads and self‐potential data. Water Resour. Res. 50(4), 3502-3522. https://doi.org/10.1002/2013WR014645. Soueid Ahmed, A., Revil, A., Steck, B., Vergniault, C., Jardani, A., Vinceslas, G., 2019. Self-potential signals associated with localized leaks in embankment dams and dikes. Eng. Geol. 253, 229-239. https://doi.org/10.1016/j.enggeo.2019.03.019. Uglebjerg, T.B., 2013. Assessment of groundwater recharge using stable isotopes and chloride as tracers (Master thesis). University of Copenhagen, Copenhagen, Denmark. van Genuchten, M. T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J, 44(5), 892-898. https://doi.org/10.2136/sssaj1980.03615995004400050002x. Vásquez, V., Thomsen, A., Iversen, B.V., Jensen, R., Ringgaard, R., Schelde, K., 2015. Integrating lysimeter drainage and eddy covariance flux measurements in a groundwater recharge model. Hydrol. Sci. J., 60(9), 1520-1537. https://doi.org/10.1080/02626667.2014.904964. Voytek, E. B., Barnard, H. R., Jougnot, D., Singha, K., 2019. Transpiration‐ and precipitation‐induced subsurface water flowobserved using the self‐potential method. Hydrol. Process. 33(13), 1784-1801. https://doi.org/10.1002/hyp.13453. Vrugt, J. A., 2016. Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation. Environ. Modell. Softw. 75: 273-316. https://doi.org/10.1016/j.envsoft.2015.08.013. Waxman, M. H., Smits, L. J. M., 1968. Electrical Conductivities in Oil-Bearing Shaly Sands. Soc. Pet. Eng. J., 8(2), 107-122. https://doi.org/10.2118/1863-A.56