Self-potential investigations of a gravel bar in a restored river corridor

. Self-potentials (SP) are sensitive to water ﬂuxes and concentration gradients in both saturated and unsaturated geological media, but quantitative interpretations of SP ﬁeld data may often be hindered by the superposition of different source contributions and time-varying electrode potentials. Self-potential mapping and close to two months of SP monitoring on a gravel bar were performed to investigate the origins of SP signals at a restored river section of the Thur River in northeastern Switzerland. The SP mapping and subsequent inversion of the data indicate that the SP sources are mainly located in the upper few meters in regions of soil cover rather than bare gravel. Wavelet analyses of the time-series indicate a strong, but non-linear inﬂuence of water table and water content variations, as well as rainfall intensity on the recorded SP signals. Modeling of the SP response with respect to an increase in the water table elevation and precipitation indicate that the distribution of soil properties in the vadose zone has a very strong inﬂuence. We conclude that the observed SP responses on the gravel bar are more complicated than previously proposed semi-empiric relationships between SP signals and hydraulic head or the thickness of the vadose zone. We suggest that future SP monitoring in restored river corridors should either focus on quantifying vadose zone processes by installing vertical proﬁles of closely spaced SP electrodes or by installing the electrodes within the river to avoid signals arising from vadose zone processes and time-varying electrochemical conditions in the vicinity of the electrodes.


Introduction
The self-potential (SP) method is a passive geophysical method, in which natural spatial and temporal variations in the electrical potential field are measured on the surface of the earth or in boreholes. The resulting SP maps and monitoring data are sensitive to flow processes in the subsurface (e.g., Doussan et al., 2002;Rizzo et al., 2004;Suski et al., 2006), but interpretation of field measurements is challenging. First, many different phenomena (e.g., water fluxes in the vadose or saturated zone, gradients in chemical potential, or redox processes) can create SP signals and it is often unclear which source types will dominate the response at a given site. Second, accurate modeling of SP responses (to given source currents) can only be achieved when detailed knowledge about the electrical conductivity distribution is available (see Chapter 4 in Minsley, 2007). Third, the selfpotential method is a potential field method and the inverse problem of retrieving the source-current distribution in the subsurface is plagued by non-uniqueness. Even if the electrical conductivity distribution of the subsurface is known, there exist an infinite number of source current distributions that can explain the data equally well. This leads to a situation in which the SP method is useful on a case-by-case basis and that its applicability at a certain site can often only be reliably assessed after the data have been acquired. Interpretation is further complicated by electrode responses that are affected by temperature variations (e.g., 0.22 mV K −1 for the electrodes used in this study; Petiau, 2000) and non-linear drift terms that are related to electrode design and age, as well as changing electrochemical conditions in the vicinity of the electrodes.
Despite the complications mentioned above, the SP method continues to receive considerable interest in hydrogeology as SP data are sensitive to contaminant transport (e.g., Maineult et al., 2004Maineult et al., , 2005Revil et al., 2009), redox processes (e.g., Linde and Revil, 2007), flow in saturated (e.g., Maineult et al., 2008;Bolève et al., 2009;Jardani et al., 2007) and unsaturated porous media (e.g., Thony et al., 1997;Doussan et al., 2002;Linde et al., 2007a), flow in fractures (e.g., Wishart et al., 2006), the water table elevation (e.g., Fournier, 1989;Revil et al., 2003;Rizzo et al., 2004), or the thickness of the vadose zone (Aubert and Yéné Atangana, 1996), etc. For example, the SP method might potentially be used to estimate water fluxes in the vadose zone (Thony et al., 1997). Such fluxes are difficult to measure in the field and are commonly inferred indirectly by differencing water content measurements over time (e.g., Vereecken et al., 2008). This richness of application areas mirrors the main limitation of the SP method; many different processes contribute to the measured response.
Self-potential source generation and the modeling of the resulting SP field are well understood under saturated conditions (Sill, 1983;Sheffer and Oldenburg, 2007). A better understanding for multiphase conditions has developed in recent years through both theoretical and experimental work (e.g., Linde et al., 2007a;Revil et al., 2007;Jackson, 2010;Allègre et al., 2010), but there is still room for improvements. These recent findings suggest that relationships between unsaturated water flux or the local hydraulic pressure gradient and SP gradients are more complex than suggested by Thony et al. (1997) and Darnet and Marquis (2004).
To decrease non-uniqueness in the interpretation of SP data, a promising approach is to treat the data as being dependent on the state of a model that describes the variations of the hydrogeological and geochemical variables of interest. This adds complexity to the problem, but makes it possible to use SP data together with other data (e.g., timeseries of hydrological head and tracer concentrations) to constrain hydrological boundary conditions, hydraulic conductivity structure, or vadose zone flow properties within a hydrogeological inverse modeling framework. This is only possible if the source contributions of different hydrological processes are accurately modeled and understood, which highlights the importance of having access to different types of hydrological and geophysical data.
A good understanding of physical processes at a given site, together with advanced signal processing and modeling, appears thus to be the only way to reliably assess the origins of the dominant contributions to the measured signals. To assess the sensitivity of SP data to hydrological processes at a restored river corridor, we performed SP mapping and monitoring on a gravel bar within a restored reach of the Thur River in northeastern Switzerland. SP monitoring in river environments is not new, for example, it has been used to study fluctuations in the water table in the vicinity of Columbia River, Washington (Timothy Johnson, personal communi- cate, 2010) and to optimize pumping rates for bank filtration at the Russian River, California (Gasperikova et al., 2008), but no conclusive results have been published to date.
Our first results concern the SP mapping and subsequent inversion of these data. We then perform an exploratory analysis of SP and hydrological time-series using wavelet analysis. For suitable periods, we compare SP time-series with those of water table position, precipitation, moisture content, and temperature. We then model, for a simplified geological model, the expected response related to the two processes that are the most likely source of the observed SP variability, namely (1) fluctuations in the water table and (2) infiltration following precipitation events.
Understanding the origins of SP signals at this site is important as one could make inferences about the hydraulic diffusivity of the aquifer if the water table effect dominates or one could obtain information about soil properties if infiltration processes dominate. SP data could then be used to evaluate, by comparing the data with those acquired at neighboring unrestored sites, the effect of river restoration on hydrological subsurface processes. Another more general motivation is that better understanding of near-surface SP sources can help to remove SP signals of shallow origin when investigating deeper phenomena, such as, volcanic activity (Friedel et al., 2004), earthquake precursors (Park, 1983), or processing magnetotelluric data (Perrier and Morat, 2000).

Thur River field site
The Thur River in northeastern Switzerland (see Fig. 1) is the largest Swiss river without natural or artificial reservoirs. It is a peri-alpine tributary of the River Rhine with a catchment area of about 1750 km 2 . River hydrology shows a typical nivo-pluvial regime. Water level and discharge variations in the Thur are characteristic of unregulated alpine rivers where neither lakes nor reservoirs attenuate the discharge. The Thur aquifer consists mainly of glacio-fluvial sandy gravels (5-7 m thick in our study area; average hydraulic conductivity 5 × 10 −3 m s −1 inferred from pumping tests; Baumann et al., 2009) overlaying thick lacustrine clays that can be considered impervious. The top of the aquifer is formed by fine alluvial sediments (fine-sand with silt). Like many other rivers, the meandering Thur was channelized at the end of the 19th century for flood protection purposes and to gain arable land (Lacey, 1930;Brookes, 1988). In an attempt to combine flood protection with ecological objectives, a more natural river environment was restored at a 5 km long reach of the Thur at Neunforn, starting in 2002. The effects of these restoration efforts are currently being investigated within the RECORD project (for details see http://www.cces.ethz.ch/ projects/nature/Record; Schneider et al., 2011).
While the channelized river was practically flowing straight prior to restoration, the riverbed morphology has changed dramatically since the removal of the northern bank stabilization and overbanks in 2002 (Trush et al., 2000;Soar and Thorne, 2001;Pasquale et al., 2010). This large widening forced the river to deposit its sediments in a typical alternate bar pattern (Tubino and Seminara, 1990). By 2005, one of those gravel bars (see Fig. 2a) had developed on the northern shore of the river with a surface exposure that depends strongly on the varying river discharge Q. The vegetation cover in Fig. 2a indicates topographic highs that reach 1.5 m above the river level at low flow conditions. At such conditions (Q ∼ 20 m 3 s −1 ) a low-lying part with bare gravel at the surface is exposed, at intermediate flow conditions (Q ∼ 100 m 3 s −1 ; 0.8 m increase in river stage compared with low flow) this part is flooded and the exposed surface of the gravel bar consists mainly of gravel covered by up to 1 m thick recent deposits of sandy loam colonized mainly by canary reed grass (Phalaris arundinacea), while the entire gravel bar is flooded at high-flow conditions (Q > 200 m 3 s −1 ; 1.4 m increase in river stage compared with low flow). This gravel bar is the focus of the SP experiments presented here.
Ten piezometers instrumented with loggers (temperature, electrical conductivity, and pressure with a sampling rate of 15 min) were installed on the gravel bar to investigate rivergroundwater interactions (Schneider et al., 2011;Vogt et al., 2010a, b). Sensors for water content (Decagon EC-5 and EC-TM) and temperature (Decagon EC-TM) were installed at different locations at 0.1, 0.5 and 1 m depth in the soil. Data were acquired with a sample rate of 30 min using a Decagon EM50 data logger (Decagon Devices Inc., Pullman, WA, USA). Raw values of soil dielectric permittivity were converted to volumetric water content using specific calibration curves for the soils at the site. A complete meteorological station (Campbell Scientific) was installed to monitor micrometeorological variables (Pasquale et al., 2010). The meteorological station includes two sensors for air temperature and relative humidity with a 10 min sampling rate (at 2.5 m and 8 m above soil level), a complete solar radiation device, two wind flow meters and an atmospheric pressure sensor. A pluviometer (OTT) with a sampling rate of 1 minute completes the station. Ground penetrating radar (GPR) and electrical resistivity tomography (ERT) geophysical data have also been acquired on the gravel bar to delineate the subsurface aquifer structure (Doetsch et al., 2011;Schneider et al., 2011).

Governing SP equations
The total electrical current density j (Am −2 ) is given by (Sill, 1983) where σ (Sm −1 ) is the electrical conductivity, E (Vm −1 ) is the electric field E = −∇ϕ, ϕ (V ) is the electrical potential, and j S (A m −2 ) is the source current density. Equation (1) is a generalized Ohm's law and Eq.
(2) is the conservation equation in the low-frequency limit of Maxwell's equations. These two equations can be combined for any electrical conductivity distribution and boundary conditions to solve for the distribution of ϕ given knowledge of j S . The total source current density for the three dominating SP sources in hydrological applications can be described by Arora et al., 2007) whereQ v (C m −3 ) is the effective charge per unit pore volume that can be dragged by the flow of the pore water; u is the Darcy velocity (m s −1 ), k b (1.381 × 10 −23 J K −1 ) is the Boltzmann constant, T (K) is the temperature, q i (C) is the charge of ionic species i dissolved in the water, t i (−) is the corresponding microscopic Hittorf number (i.e., the fraction of electrical current carried by species i in the water phase), {i} is the corresponding activity, and E h (V) is the redox potential. The first contribution is associated with the drag of excess charge in the diffuse Gouy-Chapman layer within the electrical double layer caused by the movement of the pore water. This contribution forms the streaming current and the formulation is valid under both saturated and unsaturated conditions. The second contribution is related to chemical gradients in the pore water that gives rise to diffusion currents. The third contribution is related to redox processes that only occur for the rare conditions when a path of electronic conduction (e.g., an ore body, a metallic pipe) links parts of the subsurface with different redox potentials. Revil and Leroy (2004) relateQ v at saturation to the voltage coupling coefficient at saturation C sat (V Pa −1 ) through where K is the hydraulic conductivity (m s −1 ), ρ w is the density of water (kg m −3 ), g is the acceleration of gravity (m s −2 ), and σ sat (S m −1 ) is the electrical conductivity of the saturated porous media. Experimental data suggests that C sat is largely controlled by the electrical conductivity of the pore water σ w (Revil et al., 2003) and laboratory measurements of aquifer material are rather straight-forward (Suski et al., 2006). Jardani et al. (2007) present experimental data that suggest thatQ v decreases with increasing hydraulic conductivity. Linde et al. (2007a) suggested thatQ v under unsaturated conditions should, to a first-order, scale inversely with the water saturation S w as This parameterization is similar to the one of Waxman and Smits's (1968) for a closely related parameter. This relationship is also discussed in Revil et al. (2007) and has been used by Jougnot et al. (2010). SP measurements are performed with respect to the electrical potential of a reference electrode. The measured SP data ϕ meas i (t) of electrode i at time t is given by where ϕ i (t) and ϕ ref (t) are SP responses related to the hydraulic and/or geochemical forcing defined in Eqs. (1-3). The measurements will also be affected by temporal variations in electrode coupling, electrode temperature, electrode age, and geochemistry and geology in the intermediate vicinity of each electrode. These behaviors are described by ϕ elec i (t) and ϕ elec ref (t). Despite important improvements in electrode design (Petiau, 2000), these drift terms often contaminate the long-term behavior of SP monitoring signals. Common practice is to make a linear drift correction based on SP measurements performed with the two electrodes in contact with each other in the beginning and end of the monitoring period.

SP source current inversion
The classical SP inverse problem consists of determining the position and magnitude of SP sources that can explain the observed data within the measurement errors for a given electrical conductivity distribution, while honoring any constraints on the source distribution. The SP field can be calculated by inserting Eq. (1) in Eq.
(2) to get where s (Am −3 ) is a source distribution term given by ∇ · j S ; see Eq.
(3). For a given conductivity distribution and boundary conditions, Eq. (7) can be expressed as where K is a sparse linear matrix operator that contains all information about the dimensions of the model, the electrical conductivity distribution, and boundary conditions, while ϕ is the electrical potential and s the source term at all model cells.
The inversion recovers a source model that minimizes the error between the measured ϕ obs and predicted SP response in a least-squares sense. The objective function consists of a data misfit term and a model regularization term, where P is a selector matrix that picks out the rows of K † (i.e., the inverse of K) that correspond to a SP measurement location where a potential was measured. The regularization weight λ controls the trade-off between the data misfit and the source current distribution, here quantified by the norm of the regularization operator W acting on the source model. The operator W incorporates inverse sensitivity scaling information that accounts for rapidly decaying sensitivities with distance from the measurement locations. Additionally, W promotes source solutions that are spatially sparse (Portniaguine and Zhdanov, 1999;Minsley et al., 2007); that is, it favors solutions with the fewest number of non-zero source amplitudes while still fitting the data. This sparsity constraint is non-linear and requires an iteratively re-weighted least-squares (IRLS) inversion approach. The first iteration of the IRLS inversion corresponds to the sensitivity-weighted minimum length solution, which results in models that display very smoothly varying source current distributions even if the true source locations have a limited volumetric extent. In subsequent iterations, W is updated with weights that favor models that occupy a small volume instead of being spatially smooth. We refer to Minsley et al. (2007) for more details about SP source current inversion. Even if not considered here, it is relatively straightforward to incorporate constraints from hydrogeological models or characterization concerning the regions in which source currents are expected to be located.

Wavelets
Exploratory analysis of SP field data should consider that there are many different processes that can create measurable SP signals at different temporal and spatial scales. A very useful and widely used approach to analyze non-stationary geophysical time-series is wavelet transforms (e.g., Kumar and Foufoula-Georgiou, 1997;Torrence and Compo, 1998;Grinsted et al., 2004;Henderson et al., 2009). In contrast to Fourier transforms that focus on the frequency content of a given time-series, wavelet transforms describe the temporally varying frequency content of a signal. They can also be used to evaluate how the frequency components of two different signals relate to each other over time. It appears thus that wavelet analyses could be useful to disentangle the relative contributions of the different source currents to the observed SP signals. Wavelet techniques have been used rather widely to determine the source current distribution from SP mapping surveys (e.g., Gibert and Pessel, 2001;Saracco et al., 2004), but not for monitoring purposes. In one of the rare wavelet applications to SP monitoring data, Friedel et al. (2004) analyzed data acquired at Merapi Volcano, Indonesia and found that many of the observed anomalies were associated with precipitation events.
A short description of wavelet theory is given below. We refer to Torrence and Compo (1998) and Kumar and Foufoula-Georgiou (1997) for a more detailed treatment and primary references. The wavelet used in this work is the commonly used complex valued Morlet wavelet where η is a non-dimensional time and ω 0 is the nondimensional frequency. We use ω 0 = 6 as it provides a good balance between time and frequency localization (Grinsted et al., 2004). The continuous wavelet transform (CWT) of a discrete sequence x n with a uniform time-sampling δt is defined as the convolution of x n with a scaled and translated version of ψ 0 (η) as where ( * ) indicates the complex conjugate. By varying the wavelet scale s and translating along the localized time index n it is possible to construct an image of the amplitude as a function of scale (when s increases the wavelet becomes more spread out and takes only long-term behavior of x n into account) and of how this amplitude varies with time.
The wavelet in Eq. (11) is normalized at each scale to have unit energy (Torrence and Compo, 1998). The convolution of Eq. (11) is in practice solved using Fast Fourier Transforms. The CWT has edge artifacts and it is therefore useful to define a cone of influence that identifies the regions that are sensitive to such artifacts. To assess the reliability of any features in the CWT one needs to perform a statistical test with respect to the CWT of a random process, typically red noise that can be modeled with a first order regressive process. It is then possible to test at a 5% significance level if the null hypothesis that the observed signal constitute red noise holds. The regions in which the null hypothesis can be refuted are then highlighted in the displayed CWT.
The wavelet transform offers also the possibility to create a filtered time-series x filt n between two scales j 1 and j 2 by where C δ = 0.776 is a reconstruction factor that is specific to the choice of wavelet, ψ 0 (0) = π/4, and δj is the scale sampling (12 samples for each decade is used in this study).
As W X n (s) is complex valued, it is easier to display the real-valued wavelet power spectrum W X n (s) 2 . To compare wavelet transforms of different time-series x n and y n one can calculate the cross-wavelet spectrum (WPS) as W XY n (s) = W X n (s)W Y * n (s), where W Y * n (s) is the complex conjugate of W Y n (s) and W XY n (s) is the corresponding cross-wavelet power. The wavelet coherence is defined as which can be seen as a localized correlation coefficient that varies between 0 and 1 in time-frequency space. Equation (13) is always 1 (see definition of W XY n (s)), which can be avoided when calculating the coherency by first smoothing the different contributions in time and space according to Torrence and Webster (1998). This modified coherence is the preferred measure for significance testing compared with the cross-wavelet power, which can display high values due to changes in one of the time-series only (Maraun and Kurths, 2004 Fig. 3. SP monitoring area with SP monitoring (SP1-15) and reference (SPref) electrodes, piezometers, and soil monitoring stations. The tower on which rainfall intensity was measured is located approximately 100 m away from this area.

Self-potential mapping
The SP mapping survey was carried out on 7 March 2008 using Pb-PbCl 2 -NaCl electrodes, so-called Petiau electrodes (Petiau, 2000; PMS9000 from SDEC) and a 40 MOhm impedance voltmeter). Figure 2a indicates the area that was surveyed. Measurement profiles were oriented perpendicular to the river shore with measurements every 3 m and a profile spacing of approximately 5 m (see dots in Fig. 2b). Five measurements were acquired in the vicinity of each measurement point (each reading was made after some tenths of seconds to allow stable measurements) and the median value was chosen for later processing. The measurements were performed at approximately 5 cm depth. The total data set consisted of 246 measurement points with respect to a reference electrode located in a loamy ditch. This position was chosen to obtain good and stable electrical coupling conditions. The data were corrected for a linear drift of 2.8 mV over time using five drift measurements (i.e., the SP signal is recorded when the measurement electrode is located close to the reference electrode) acquired during the day. For further interpretation, the reference (0 V) was assigned to the position indicated in Fig. 2b. The data were detrended using linear regression in the direction of water flow and interpolated using ordinary kriging (Deutsch and Journel, 1998). The kriging was performed with a spherical model that fitted the detrended experimental semi-variogram with a nugget of 3 (mV) 2 , an effective range of 18 m and a sill of 20 (mV) 2 (see Linde et al. (2007b) for a more detailed description of kriging of SP data). Figure 2b displays the kriged map together with the linear trend model. The vegetated region with sandy loam soils shown in Fig. 2a have values in the range of −8 to −13 mV, whereas the regions in the lower-lying bare gravel have values close to 0 mV. One possibility is that the source regions are mainly located within the finer and thicker soils instead of the exposed gravel, though another possibility is that the SP sources originate from groundwater flow and aquifer heterogeneity in the gravel aquifer itself that happen to coincide with the soil boundaries seen on the surface. However, the strong SP gradients in Fig. 2b indicate that the sources are at least partly located in the shallow subsurface.
A repeat survey (not shown here) was performed on 4 March 2009, which was carried out during a day when the river and groundwater level was 45 cm higher than during the previous survey. The resulting magnitudes are close to five times lower, which strongly suggest that the origin of the SP signals are more related to vadose zone processes than groundwater flow.

SP inversion
We use the 3-D SP source current inversion method of Minsley (2007)  The discretization was 1 m in the x-and y-direction, whereas it was 0.5 and 1.0 m in the vertical direction for the gravel and clay layers, respectively. The inversion was first carried out until the model reached a mean data misfit of 1.3 mV before an additional iteration was carried out with compactness constraints. Figure 2c and d display the magnitudes of the source currents between 1.5-2 m and 3.5-4 m depth, respectively. They clearly indicate zones of negative source currents that approximately cover the outer part of the soil-covered region of the gravel bar. The maximum amplitudes of the source currents are found at 3 m depth. The magnitude of the source current decreases with deeper depths and is negligible at 5.5-6 m depth (Fig. 2e). The inversion results thus indicate that the source currents are found in the shallow subsurface at locations that mainly correspond to the soil-covered part of the gravel bar. Note that the SP inversion cannot resolve the accurate depth of the sources due to the inherent nonuniqueness of the inverse problem and the simplicity of the resistivity model.

SP time-series and wavelet analysis
The self-potential monitoring was performed from 13 February until 11 April 2009 in the region indicated in Fig. 2a. The locations of the SP monitoring electrodes, piezometers, and soil monitoring sites (temperature and water content) are shown in Fig. 3. A total of 16 Petiau electrodes were installed at 20 cm depth from the surface. A thin layer of fine sediments was added to the electrode-soil interface to improve the electrical contact and the region around the electrode was wetted before the monitoring started. The reference electrode was located in a thick loam layer. We used a Campbell Scientific Inc. CR1000 data logger that measured and stored the voltage difference between the measuring electrodes and the reference electrode every 5 s. A 12 V battery was used to power the CR1000. The connections between the electrodes and the data logger were achieved using isolated, but unshielded, copper wire. The data were processed by median filtering (no detrending) and sampled every 15 min for subsequent analysis. Figure 4a shows time-series of rainfall intensity during the monitoring period, in which two periods of more significant rainfall occurred between days 60-70 (first rainy period) and 82-87 (second rainy period). Figure 4b shows time-series of water content at 10 cm depth (see Fig. 3) with increases in water content corresponding to the rainfall periods. The hydraulic head data of Fig. 4c have a delayed response to the rainfall and display some uncorrelated events that are attributed to snowmelt and rainfall upstream. The soil temperature data at 10 cm depth (Fig. 4d) show that the ground was partly frozen until day 57 followed by a successive warming of the soil with daily fluctuations of a few degrees. Selected SP data (Fig. 4e) show that the main periods of SP variability correspond to the periods of rainfall. One can also observe substantial long-term variations in SP11, which is most likely related to electrode drift. No SP data were acquired in the beginning of the second rainy period as interfering ERT geophysical measurements were acquired during this period.
The WPS of the precipitation data ( Fig. 5a) display significant energy covering the shortest period of 30 min to a period of approximately two days. The WPS of the water content (Fig. 5b) shows a similar pattern, but with less high-frequency content and limited response to the rainfall between days 45 and 55, possibly due to frozen ground conditions. The WPS of the hydraulic head (Fig. 5c) displays energy during the two main periods of rainfall, but there is also significant energy with periods longer than a day corresponding to variations in the hydrological conditions in the upstream region of the catchment. The WPS of the soil temperature data (Fig. 5d) shows only well-defined variability associated with daily fluctuations. Four different SP electrodes were chosen with WPSs that visually appear to be the most related to the state variables discussed above. The WPSs in rainy periods, whereas (h) SP2 has a different behavior and shows a rather significant daily variation, which is attributed to the fact that SP2 is largely insensitive to what happens during the rainy periods. The wavelet coherencies in Fig. 5 were calculated for the (i) rainfall intensity and SP5, (j) the water content and SP3, (k) hydraulic head and SP11, and (l) temperature and SP2 following Grinsted et al. (2004). The coherent events indicate in most cases time-varying phase relations. The more significant coherent periods with consistent phase relations occur for the precipitation and water content data for periods of 1 to 4 days during the two more rainy periods.
A close up of the second rainy period is shown in Fig. 6. The long-period contributions of more than two weeks have been removed using Eq. (12) to focus on short-term variability without superimposed long-term trends or drifts. The rainfall intensity, detrended water content, and detrended hydraulic head variations are shown in Fig. 6a-c. Detrended SP time-series are shown in Fig. 6d-i. It appears that (d) SP11 and (g) SP6 are very sensitive to the rainfall intensity, such that even short periods of rainfall within the rainy periods have clearly defined SP peaks. The variability in (e) SP3 shows a rather close resemblance with variations in water content, whereas (h) SP1 show more of a gradual buildup of the SP signal over time with an abrupt decrease after the end of the rainy period. Finally, both (f) SP8 and (i) SP12 appear to be mainly correlated with the hydraulic head. The correlations observed in Fig. 6 are similar to those obtained in the first rainy period (not shown). These examples indicate that the dominating SP source-generating processes can vary dramatically within a rather small monitoring area. We found no significant relationship between the degree of correlation between the SP signals and water table fluctuations as a function of the thickness of the vadoze zone (not shown here), indicating that the SP signals are likely more related to soil heterogeneity than water table dynamics.

SP modeling
An idealized geological model was created to investigate how SP signals are expected to vary with respect to a rising water table and precipitation (Fig. 7). Two simulations of each type are investigated. Tests 1 and 2 simulate a water table increase with a rate of 10 cm h −1 during one hour from an initial level H0 of 0.5 and 1.5 m depth, respectively, which represents a typical rate in response to a moderate precipitation event upstream. Tests 3 and 4 simulate two rainfall events with an intensity of 0.3 mm h −1 during 1 h with a constant water table at 0.5 m and 1.5 m depth, respectively. The simulation time is 2 h for all tests, where the hydrological events take place during the first hour and the relaxation of the SP signal is investigated in the second hour (Fig. 7b, c and d).   The model geometry (Fig. 7a) consists of an upper 1 m thick layer with gravel on the left side in which we investigate variations of the SP signal, and loam on the right side, in which we have placed our reference (i.e., the electrical potential is zero). The reference is located in the loam because it is the most stable location over time and because it corre-sponds to the field situation. A uniform gravel aquifer is located between 1 and 9 m depth followed by a 20 m thick clay aquitard. The boundary conditions for the electrical problem is electrical insulation at all boundaries. Hydrological boundary conditions are no flow boundaries at the sides and imposed pressure at the clay-gravel aquifer interface for all simulations (Fig. 7c-d). A no flow boundary at the top of the gravel is defined for the increasing water table simulations, while a prescribed flux into the gravel is prescribed for the precipitation experiment as displayed in Fig. 7b. No flow boundaries are imposed on top of the loam as it is assumed that infiltration is negligible in the loam for the examples and time scales considered here. The choice of Neumann boundary conditions for the electrical and hydrogeological problems on the sides can be motivated by symmetry arguments. The aim of the modeling is to investigate the perturbation of the SP field caused by the soil heterogeneity, while modeling an essentially infinitely sized aquifer system.
The modeling is performed in 2-D using finite element calculations in COMSOL Multiphysics 3.5 with a mesh consisting of close to 50 000 triangular elements with specific mesh refinement in the unsaturated zone. The electrical resistivity of the clay is 25 m, whereas those of the loam and gravel are modeled by with the electrical formation factor F being 12 for the gravel and 4 for the loam, Archie's saturation exponent n = 2, the initial electrical conductivity of the groundwater was σ w = 0.04 S m −1 for both soils, while the electrical conductivity of the rainwater is 0.002 S m −1 . A surface conductivity σ s = 0.002 S m −1 was chosen for the loam, whereas it was assumed to be zero for the gravel. We used the parameterization of van Genuchten (1980) for the relative hydraulic conductivity K r and capillary pressure P c functions as where S w is the water saturation, S e is the effective and S wr is the residual water saturation, respectively, and m and α (m −1 ) are soil-specific parameters. We used typical parameters for gravel and loam (e.g., Carsel and Parrish, 1988) as outlined in Table 1. The hydrological problem is solved using Richard's equation. Laboratory measurements using sediments retrieved from neighboring cores suggest that C sat is −21 ± 3 mV m −1 for measurements performed at 20 • C and with σ w = 0.034 S m −1 . However, these cores were severely disturbed and it was impossible to obtain representative estimates of σ s and K, which also indicate that our C sat estimates might be biased. For the modeling, we decided instead to assume thatQ v (S w = 1) is 0.48 C m −3 for the gravel and 28 C m −3 for the loam following the experimental relationship of Jardani et al. (2007) for a saturated medium. We use Eq. (5) to modelQ v (S w ).
The SP signal for test 1 (Fig. 8a) displays a positive SP signal over time that is at its maximum after 1 h and then decreases slowly as a certain upward flow continues to occur in the vadose zone to reach hydrostatic equilibrium (Fig. 8b). At a given time, the SP signal is constant with depth throughout the part of the vadose zone in which no changes in water content occurs. No SP signals (Fig. 8c) occur in this region for test 2 for which the source currents at the water table are the same throughout the model. This effect occurs as the water table (Fig. 8d) is located in a region of uniform geological media. These results indicate that SP signals in the unsaturated zone only occur when there are lateral contrasts in the sediments in which the water table rise or where the water content changes with time. It appears thus that any empirical relationship between water table dynamics and SP signals will vary between electrode positions and that it might be highly nonlinear (e.g., no sensitivity at all until a contrast is reached). Heterogeneity will thus play a key role in determining the SP signals associated with water table fluctuations both in terms of determining the location and magnitude of source currents, but also in determining the electrical conductivity distribution that also strongly affects the SP magnitudes; see Eq. (7).
The SP response (Fig. 9a) to rainfall in test 3 indicates that infiltration creates a vertical SP gradient in the vadose zone with negative magnitudes within the vadose zone (Fig. 9b). A similar SP response (Fig. 9c) is shown for test 4, but with lower magnitudes. The magnitudes of the SP signals are much lower than those observed experimentally by Doussan et al. (2002), which can partly be attributed to different soil properties. Even if ignored here, there might be significant contributions from diffusion currents -see Eq.
(3) -as the infiltrating rainwater has a much lower ionic content.

Discussion
Self-potential mapping makes it possible to investigate the overall SP field on the surface at high spatial resolution, but the interpretation needs to consider that contributions of SP sources located in the vadose zone might vary over short spatial and time scales. We suggest that SP mapping for hydrogeological purposes should be performed following a relatively long period (e.g., a week) of stable meteorological and hydrological conditions to minimize the effects of sources in the vadose zone. SP monitoring can be useful to monitor vadose zone processes and to better understand the interrelations between different processes and mechanisms of SP source generation. The monitoring is however affected by long-term electrode drifts and what appears to be a high sensitivity to electrochemical conditions in the close vicinity of the electrodes and possibly wetting conditions at the electrode tip. Repeated (e.g., weekly) SP mapping with measurements in the vicinity of the monitoring electrodes might be useful to remove some of the problems with monitoring data. Even so, the data acquired in this study raise doubts about the usefulness of slowly varying (e.g., periods of weeks) natural variations in SP signals when their amplitudes are in the range of some 10 mV and electrodes are placed in the vadose zone.
That the SP method is cheap, light, and fast is often presented as making it particularly useful in regions where accessibility is limited and where limited data are available. We argue that for anything but qualitative applications it is crucial to have access to supplementary data about the geology, depth to water table, water chemistry, and meteorological conditions. It is also necessary to have a reliable model of the electrical conductivity distribution over time.
The SP time-series at the Thur River site display a non-stationary behavior with varying sensitivities to possible forcing parameters (rainfall intensity, water table, etc.). Time-series analysis with wavelets is therefore a very suitable tool to better understand casual and sometimes intermittent relationships between SP signals and other state variables.
Our results suggest that most of the SP signals on the surface of the gravel bar are related to sources in the vadose zone, which prohibit attempts to use SP signals to infer flow patterns in the aquifer. The magnitudes of the observed signals are much higher than those obtained by modeling, which is likely due to that the scaling relationQ v (S w ) suggested in Eq. (5) is too simple. Note that approaches that implicitly assume thatQ v (S w ) = 1 (Perrier and Morat, 2000;Darnet and Marquis, 2004) would predict even smaller SP magnitudes for the same choice of material properties. The high sensitivity to the heterogeneity of vadose zone states and fluxes forms in our mind an important motivation for continued SP research. In fact, measurable SP signals originating in the vadoze zone will only occur between electrodes located at the same depth in cases of differences in vertical water fluxes, lateral fluxes, or heterogeneity in electrical and soil properties. To advance understanding, we suggest that vertical profiles of SP electrodes should be installed in wellcharacterized and well-instrumented soils, such that more realistic modeling can be performed than what is presented here. Such time-series could allow us to understand how to infer fluxes in the vadose zone from SP measurements, which would necessitate an accurate soil-specific function to predict Q v (S w ).
Another potentially fruitful approach would be to perform SP monitoring and mapping within river and lake systems (i.e., under saturated conditions) to investigate infiltration and exfiltration processes without having to deal with the highly complex vadose zone response while assuring good electrical coupling conditions. This would necessitate having access to not only river temperature and electrical conductivity information, but also that installations are made at locations at which deposition and erosion processes are negligible on the time-scale of the monitoring period.
It might appear surprising that the groundwater flow component to the SP signals is so low at our study site, but the reason for this is very simple. Groundwater flow is taking place in a rather thin resistive aquifer (5-7 m thick), while insignificant groundwater flow takes place in the thick electrically conductive underlying alluvial clay. SP theory (e.g., Sill, 1983) dictates where source currents are located within the gravel. However, the resulting SP distribution is not only dependent on the source current distribution, but also on the electrical conductivity distribution within and outside of the source region. In fact, the alluvial clay channels the electrical current, which decreases the resulting SP signal drastically compared to the case of a resistive basement. This comes naturally from the boundary condition of Maxwell's equations that state that the tangential component of the SP gradient (i.e., minus the electrical field) is equal on both sides of a boundary with different electrical properties (e.g., clay and gravel in this case). The SP signals related to groundwater flow would be about one order of magnitude more important if the aquifer was underlain by a resistive bedrock.

Conclusions
We find that wavelet transforms can be extremely useful to disentangle the temporally varying interrelations between SP monitoring signals and other state variables (e.g., rainfall intensity, soil temperature, hydraulic head, water content variations) and thus better understand the main SP source mechanisms at a given site. At the Thur River it appears that the SP signals are mainly determined by local soil heterogeneity together with variations in water content, infiltration, and groundwater level. Contributions from groundwater flow appear to be of limited importance mainly due to the conductive underlying alluvial clay.
We suggest that future SP monitoring experiments in river environments should focus on either (1) estimating water fluxes in dedicated and well-instrumented soil profiles (several SP electrodes installed at different depths in addition to temperature and soil moisture sensors) or (2) mapping and monitoring on the river bed to avoid the influence of the spatially and temporally varying non-linear SP signals associated with vadose zone processes.
Approaches that interpret SP signals only in terms of natural water table fluctuations can only be expected to work in environments with no or very limited infiltration (i.e., dry regions, urbanized regions, or by installing plastic liners on the surface between the electrodes). Even if the understanding of self-potential signals in unsaturated media has increased in the last years, it appears that more dedicated field experiments and laboratory work with appropriate complimentary data are needed before SP can become a quantitative monitoring tool of hydrological processes in river corridors.