Spectral induced polarization measurements for predicting the hydraulic conductivity in sandy aquifers

Field and laboratory spectral induced polarization (SIP) measurements are integrated to characterize the hydrogeological conditions at the Schillerslage test site in Germany. The phase images are capable of monitoring thin peat layers within the sandy aquifers. However, the field results show limitations of decreasing resolution with depth. In comparison with the field inversion results, the SIP laboratory measurements show a certain shift in SIP response due to different compaction and sorting of the samples. The SIP data are analyzed to derive an empirical relationship for predicting the hydraulic conductivity ( K). In particular, two significant but weak correlations between individual real resistivities ( ρ) and relaxation times ( τ ), based on a Debye decomposition (DD) model, with measured K are found for the upper groundwater aquifer. The maximum relaxation time (τmax) and logarithmically weighted average relaxation time (τlw) show a better relation with K values than the median valueτ50. A combined power law relation between individual ρ and τ with K is developed with an expression of A · (ρ) · (τlw) , whereA, B andC are determined using a least-squares fit between the measured and predicted K. The suggested approach with the calculated coefficients of the first aquifer is applied for the second. Results show good correlation with the measured K indicating that the derived relationship is superior to single phase angle models as Börner or Slater models.


Introduction
A key prerequisite for reliable prediction of groundwater movements is the hydraulic conductivity (K).The pumping test, grain size or coring sleeve analyses have traditionally been the standard methods used to evaluate the hydraulic parameters of subsurface material for characterizing groundwater aquifers.These methods are expensive, slow and often unavailable due to disturbance during drilling to get the borehole samples.Partly for these reasons, electrical methods are playing an increasingly important role in predicting the aquifer hydraulic's parameters (e.g., Zisser et al., 2010a;and Khalil, 2012).
The relationship between direct current (DC) resistivity ρ and K is first considered as this makes up the bulk of the pertinent literatures published over the last 30 yr (e.g., Heigold et al., 1979;Niwas and Singhal, 1985).Further, the relation between the real conductivity part (σ ) and permeability was discussed and numerically investigated by many authors (e.g., Sen et al., 1981;Huntley, 1986;Bernabe and Revil, 1995;and Frohlich et al., 1996).Empirical and semiempirical relationships have been established in terms of power law relations between various aquifer parameters and those obtained by resistivity measurements (e.g., Attwa et al., 2009).Because the hydraulic parameters depend on the porosity and the geometry of the pore space, K cannot be uniquely determined by DC resistivity alone without further assumptions (Hördt et al., 2007).
The induced polarization (IP) of non-metallic minerals is generally referred to as interface or membrane polarization (Marshall and Madden, 1959;Vinegar and Waxman, 1984).
Published by Copernicus Publications on behalf of the European Geosciences Union.
In the absence of metallic conductors, spectral induced polarization (SIP) phenomena are commonly associated with polarization effects related to an electrochemical double layer (EDL), which describes the organization of ionic charges at the interface between solid and fluid (Revil and Florsch, 2010).The EDL provides the conceptual background for the electrochemical processes considered to be responsible for a large amount of the observed SIP response (e.g., Leroy et al., 2008;and Revil and Cosenza, 2010).Accordingly, both ohmic (i.e., amplitude) and capacitive (i.e., phase) parts of interface polarization contain additional information about the textural characteristics of sedimentary rocks (e.g., Schön, 1996;Kemna, 2000;Slater, 2007;and Blaschek and Hördt, 2009).
By measuring chargeability or, equivalently, the phase shift between current and voltage, information about the pore space can be gained.Previous workers have shown that the IP mechanisms can be very sensitive to changes in lithology and pore fluid chemistry (e.g., Pelton et al., 1978;Weller et al., 2010b;and Skold et al., 2011).Environmental examples for the successful use of SIP include the detection of clay units (e.g., Hördt et al., 2009;Attwa et al., 2011;and Attwa and Günther, 2012), the detection of both organic and inorganic contaminants (e.g., Olhoeft, 1984Olhoeft, , 1985Olhoeft, , 1992;;Chen et al., 2008Chen et al., , 2012;;andAbdel Aal et al., 2009, 2010) and hydraulic conductivity estimation (e.g., Tong and Tao, 2007).Although the SIP method offers potential for subsurface structure and process characterization, in particular in hydrogeophysical and biogeophysical studies (Kemna et al., 2012), it is not as widely used as other electrical methods (e.g., DC resistivity), and its full potential has yet to be realized.
The hydrogeophysical research of the IP method has shown that SIP data can be correlated with physical properties of the pore space in non-metallic soils and rocks, such as the specific surface area (S por ) and K (e.g., Pape et al., 1987;Börner et al., 1996;and Weller et al., 2010a).There is a growing interest in the use of SIP for a wide range of environmental applications, in particular those focused on hydrogeological investigations.A clear link between hydrogeological properties and IP/SIP parameters has been empirically documented by various studies (e.g., Binley et al., 2005;and Hördt et al., 2007).Binley et al. (2010), Kruschwitz et al. (2010) and Titov et al. (2010) also pointed out that the understanding of the detailed nature and origin of such linkages still lies ahead.
Recently, new models have been developed to describe the complex conductivity of saturated and unsaturated sands and clayey materials in the low frequency domain (Revil and Florsch, 2010;Revil, 2012;Revil et al., 2012aRevil et al., , 2013a)).Revil et al. (2013a) presented a comparison between recently derived petrophysical models and a new set of petrophysical measurements including permeability, complex conductivity, and streaming potential data on saprolite core samples.They proved that the quadrature conductivity (i.e., polarization) is not controlled by the diffuse layer (the outer component of the electrical double layer) but more likely by the Stern layer (the inner layer of the electrical double layer).The surface conductivity is dominated by electromigration in the electrical diffuse layer coating the pore water-mineral interface.Revil (2013) indicated that the diffuse layer is not the main contributor to low-frequency polarization of clayey materials and also that the Maxwell-Wagner polarization is not the dominating mechanism of polarization at intermediate frequencies.Accordingly, he presented an extremely simple unified model based on the low-frequency polarization of the Stern layer.This model is proposed for all sedimentary rocks, saturated and unsaturated clayey materials.It offers the possibility to predict the conductivity and permittivity (or alternatively the quadrature conductivity) as a function of frequency, clay content, and clay mineralogy.It is also used to predict permeability in the critical frequency observed at low frequencies in the quadrature conductivity.
Laboratory SIP investigations covering a wide frequency range provide more information on the spectral behavior of conductivity amplitude and phase shift.There are no universal physically based models that describe the frequency-dependent complex conductivity response of sediments and, consequently, macroscopic representations have been adopted in the literature (Binley et al., 2005).The SIP response can be described as a superposition of relaxation processes, which might reflect grain-size distribution (Lesmes and Morgan, 2001).The induced polarization decay contains much more information than single frequency or total chargeability data that are often used in traditional IP measurements.Significant progress has been made over the last decade in the understanding of the microscopic mechanisms of IP; however, integrated mechanistic models involving different possible polarization processes at the grain/pore scale are still lacking (Kemna et al., 2012).
Various conceptual and theoretical relaxation models exist to describe the SIP mechanism within the rocks and the relaxation process after terminating the current.Different approaches have been reported for the computation of the relaxation time (τ ) distribution.One of the most popular models is defined as Debye decomposition (DD) (Nordsiek and Weller, 2008).The DD can be applied to any measured spectral data sets independently of whether the spectra exhibit constant phase, by the same performance as Cole-Cole or other behavior (Weller et al., 2010b) present.Binley et al. (2005) use the mean relaxation time (τ m ) from a Cole-Cole model of the complex resistivity to estimate the permeability of sandstone.Their measurements reveal evidence of a relationship between τ and a dominant pore throat size.Revil and Florsch (2010) presented a model indicating that the permeability depends linearly on relaxation time.Revil (2012) developed a new robust model named PO-LARIS describing the complex conductivity of shaly poorly sorted sand.He showed that the permeability can be predicted inside 1 order of magnitude using the formation factor F and the quadrature conductivity.Revil et al. (2013b) developed a simple polarization model accounting for the dependence of the in-phase and quadrature conductivities on the porosity, cation exchange capacity (or specific surface area of the material), and salinity of the pore water for simple supporting electrolytes like NaCl in isothermal conditions.They showed a linear relationship between the main relaxation time and the pore throat size for clean sand and clayey sand samples.Revil (2013) predicted the permeability for various clean sands, clayey sands and sandstones using a low-frequency relaxation time τ and the formation factor F .Accordingly, these studies indicate that there is a linear relation between the relaxation time and the pore size characteristics and the permeability of the materials.
The main objective of this paper is to use the SIP field and laboratory data to predict the K of unconsolidated sandy aquifers at the hydrogeological test site Schillerslage, Germany.A practical empirical relation will be introduced to calculate K.The format of this paper is as follows: (a) SIP field (sounding and profiling) and laboratory data will be acquired for describing the geological and hydrogeological characteristics of the subsurface Quaternary aquifers.(b) Field and laboratory SIP data are used to evaluate various published approaches to predict K values.The Debye decomposition model (Nordsiek and Weller, 2008) is applied to determine the τ from available laboratory SIP data.(c) The importance of combining electrical resistivity (ρ ) and τ for computing Kis emphasized.Finally, the calculated K from field and laboratory data using various approaches are compared with those measured from the pumping test and coring sleeve information.

Material and methods
The SIP technique is a complex multi-frequency extension of the DC resistivity.Alternating currents of various frequencies are injected and the phase shift between voltage and current signals is measured in addition to the resistance (amplitude ratio).If the porous medium consists of clay-free, uniformly sized particles, the IP decay curve would consist of a single exponent decay of relaxation time constant, τ , (Vinegar and Waxman, 1984, 1987, 1988).Revil and Florsch (2010) assumed that the Cole-Cole time constant, τ , can be related to the square of the grain diameter d as where D is the diffusion coefficient in the Stern layer, which is temperature dependent and varies ∼ 2 % per degree Kelvin.Titov et al. (2010) showed the relationships between the modal pore throat diameter (R) in µm and the characteristic mean relaxation time, τ m in s.A positive logarithmic dependence was derived: R = 16.8 ln τ m − 2.15; (2) Titov et al. (2010) suggested that the polarization of the sandstones is controlled with clay content.They also showed a good correlation between R and permeability k.
The estimated geoelectrical properties of earth materials may be represented by complex electrical resistivity (ρ * ): where ρ * , ρ and ρ are the complex resistivity, real resistivity part and the imaginary part of ρ * , respectively, and ϕ is the phase angle of ρ * that can be also written as follows: The relationship between K and σ will depends on the physical and chemical properties of mineral and pore fluids and it indicates whether electrolytic conductivity (σ el ) or interfacial conductivity (σ int ) dominates measured σ .Traditionally, the real part of resistivity (i.e., DC resistivity) has been used to estimate K values for unconsolidated aquifers.For example, positive or negative linear log-log relationships between K and ρ and σ have been explained by many authors (e.g., Heigold et al., 1997;Purvance and Andricevic, 2000;Attwa et al., 2009;and Attwa, 2012), using To apply Eq. ( 5) it would be necessary to estimate a and b based on the comparison of geoelectrical measurements with measured K values (e.g., pumping tests or grain size analyses).For sandy clayey soils the correlation of K with ρ is directly proportional on a large scale, but sometimes on a local scale the correlation is inversely proportional (Mazac et al., 1990;Attwa, 2012).
There exists a substantial volume of literatures demonstrating a power dependence of the imaginary part of conductivity (σ ) on S por (e.g., Börner et al., 1996;Slater and Lesmes, 2002;Slater et al., 2006;and Weller et al., 2010b).Börner et al. (1996) observed that the phase shift (ϕ) for most of their samples was particularly constant over a broad frequency range and suggested an empirical relation to estimate K from data for one frequency only, preferably around 1 Hz.Consequently, the basic assumption of this model is that the measured phase shift of complex conductivity is independent of measurement frequency.However, for most rocks, this assumption is not valid.The proposed model (PARIS equation, Pape et al., 1987) for a constant-phase-angle (CPA) to calculate K (in m s −1 ) is where a and c exponents are adjustable parameters, F is the formation factor and S por is the specific surface area.The parameters in Eq. ( 6), suggested by Pape et al. (1987), are valid for consolidated sedimentary rocks.Börner et al. (1996) modified Eq. ( 6) to be valid for unconsolidated rocks.Whereas Pape et al. (1987) suggested a = 0.00475 for consolidated sediments, Börner et al. (1996) used a = 1 instead.Here, we work with a = 0.00475, which was also recommended by Weller and Börner (1996) and Hördt et al. (2009) for unconsolidated sediments.The exponent c was found to be in the range 2.8 < c < 4.6 (Börner et al., 1996), depending on the method of K measurements.Börner et al. (1996) presented a modification to Archie's law to calculate the F as where σ w is the pore fluid conductivity.Börner et al. (1996) suggested to use l = 0.1 in unconsolidated sediments.Note that Eqs. ( 6) and ( 7) are not universal, but they were applied for wide hydrogeological investigations (e.g., Börner et al., 1996;andHördt et al., 2006, 2009).S por was calculated after Weller et al. (2010b) for unconsolidated sandy aquifers as where S por is in µm −1 and σ in mS m −1 .
The second approach was suggested by Slater and Lesmes (2002) using the CPA model.The Slater and Lesmes (2002) model is based on a relationship of the imaginary conductivity part with the diameter of the grain corresponding to the smallest 10 percent portion of the sample (i.e., d 10 ).Furthermore, this approach is based on the inverse relation between K and σ around 1 Hz as where σ is given in µS m −1 .m and n are adjustable parameters with m = (2.0 ± 0.3)10 −4 ; n = 1.1 ± 0.2, and K is in m s −1 .It is noticed that the exponent from Slater and Lesmes model ( 2002) is smaller than in the model of Börner et al. (1996); hereinafter these models are referred to as "Slater model" and "Börner model", respectively.It is generally accepted that τ (tau) is controlled by the diffusion process (Revil and Cosenza, 2010) in the electrochemical double layer (EDL).A power law relationship between some relaxation time constants and permeability (k) has been reported in published data sets (Binley et al., 2005;Kemna et al., 2005;Tong et al., 2006a, b;Koch et al., 2011): where a and b are formation-specific parameters.Zisser et al. (2010b) established a power law relationship between the mean relaxation time (τ m ) and S por .They proved also that the relationship between the median relaxation time (τ 50 ) and S por is stronger than that between τ m and S por .Accordingly, in presence of a weak relation between the τ and the permeability and/or K, this relation cannot be applied.Weller et al. (2010a) presented a generalized power law relation using DC resistivity, chargeability (M), and mean relaxation time (τ m ), which were derived from DD, to predict k for isotropic and anisotropic samples, The four empirical parameters a, b, c and d are determined by a multivariate regression analysis.

Geological overview and measurement strategy
The test site Schillerslage is located northeast of Hannover and shows a typical geological structure for the Quaternary sediment basin of northern Germany.Based on borehole information, two sandy aquifers separated by a fine-grained till layer can be recognized (Fig. 1, right).The upper aquifer, down to a depth of about 12 m, consists of medium to partly coarse sand and thin interbedded peat layers.The aquiclude till layer with high clay content is, in general, 12 m to 16 m deep and varies also in thickness.The lower aquifer (i.e., the second one) of 5 m thickness consists of slightly limy medium grey sand and it is overlying Cretaceous marls.
The SIP measurements were started by measuring a sounding curve using a Schlumberger configuration with AB/2 values from 1.5 to 125 m and MN/2 values of 0.5 and 5 m.Additionally, a short SIP profile (P1) was measured crossing the borehole ENG_03 and the center of sounding curve (see Fig. 1).The IP profile P1 was measured using 21 electrodes and 2 m electrode spacing.The P1 profile was followed by measuring a long profile (P2) by using 36 remote units to characterize lateral heterogeneity.Moreover, electrode spacing of 7 m was used to acquire the 2-D data across the P2 profile (Fig. 1, right).
SIP laboratory measurements were done on 33 core samples at different depths (from 0.65 to 23.15 m) of wells ENG_03 and ENG_08 (close to ENG_03).Then, direct K measurements were carried out by Sass (2010) using coring sleeves (16 samples) and grain size analysis (13 samples), for the upper and lower aquifers, respectively.The coring sleeves were carried out on whole cores (8 cm diameter and 100 cm length) and soil samples (6 cm diameter and 4 cm length).For grain size analysis, we used equations after Kozeny (1927) and Beyer (1964) to derive K and found that the differences are insignificant.Both field and laboratory SIP measurements were conducted to calculate K values using various approaches.Nine laboratory SIP samples located at the same depth of those used for coring sleeves were used to derive an empirical relationship for predicting K values.The inverted field (sounding and profiling) and measured laboratory SIP data were correlated and evaluated to assess the differences in K values.
The SIP measurements (resistivity amplitude and phase angle) are conducted over a wide range of frequencies.To minimize induction effects, a dipole-dipole configuration  was used to acquire the profile's SIP data.The field and laboratory SIP measurements were carried out using SIP256C and SIP Fuchs equipment, respectively, by Radic Research (Radic et al., 1998;Radic, 2004).The four-electrode device SIP Fuchs measures the complex resistivity over almost seven decades of frequency (from 0.0457 to 12 kHz).On the other hand, the multi-electrode instrument SIP256C measures the complex resistivity over more than three decades (from 0.625 to 1 kHz).In both instruments, the remote units register the voltage or current measurements, digitize the data and transfer them to controlling PC through a fiber optic cable to avoid interference of the current supply cables with the voltage measurement.Acquisition of very low-frequency data is limited by the data acquisition time.For one profile the total measuring time was about 2 h, which can easily become excessive if even lower frequencies (< 0.625 Hz) are required, even though the system efficiently measures voltages at all receiver channels simultaneously.

Inversion results and interpretation
A Gauss-Newton algorithm using a smoothness constraint with fixed regularization was chosen for sounding and profiling inversions, which were done by using DC1dinv and DC2DInvRes softwares (Günther, 2004), respectively.The 2-D inversion algorithm was implemented for a global regularization scheme using a first-order smoothness constraint (Günther, 2004).As introduced by Constable et al. (1987), we used different weights for horizontal and vertical model boundaries such that a x = λ and a z = λw z .The regularization parameter (λ) is a weight that adjusts the degree of importance of the model smoothness constraints versus the data misfit.A small value of λ (or w z ) produces a highly structured model with huge parameter contrasts, explaining data well, whereas a big value will not be able to fit the data but provides a smooth model.
While the apparent resistivity/conductivity amplitudes are almost constant over frequency, the phase shows inductive or capacitive coupling at frequencies starting at about 10 Hz.To avoid distortion by the coupling effect, only low frequency data were used in the inversion process.The 2-D processing included the rejection of data points with bad quality from the measured data, i.e., for which the stacking error was above 1 %.In addition, 3 % error was added to the stacking error to account for systematic error components in absence of reciprocal data.

Field inversion results and interpretation
Figure 2 shows the inversion results of the IP sounding at 366 mHz.A good fit between the measured and theoretical resistivity data can be observed (Fig. 2a).On the other hand, a bad comparison can be observed at AB/2 higher than 50 m (Fig. 2b, right).This can be attributed to the electromagnetic coupling, which can be noted by abrupt changes in the phase angle readings.
The inversion results of resistivity amplitudes show a good comparison with the borehole data (Fig. 2c, left).A lowresistivity layer (< 50 m) can be noticed between the two aquifers, corresponding to the till layer.At about 23 m depth, another low resistivity layer can be observed corresponding to the Cretaceous marl layer.These low resistivity layers are characterized by high phase angle values (Fig. 2d, left).
As in Olayinka and Weller (1997), the 1-D inversion results and the borehole information were used as a starting model for 2-D inversion to constrain the inversion process.Figure 3 shows the inversion results of the profile P1.Low  frequency (0.625 Hz) SIP data were used in the inversion process.For the 2-D inversion, the regularization anisotropy was set to w z = a z /a x = 0.01 in order to achieve predominantly layered structures.It is clear from the resistivity amplitude (Fig. 3a) that the resistivity decreases with depth.In comparison with the borehole data, the heterogeneities of thin layers are not well defined from the resistivity amplitude image; whereas within the upper sandy aquifer's thin layers with high phases (> 8 mrad) and sharp boundaries they can be well noticed from the conductivity phase model (Fig. 3b), that corresponds to the peat layers.At about 13 m depth, there is a sharp boundary between low (< 4 mrad) and high phase values (> 13 mrad) which corresponds to the boundary between the second sandy aquifer and the till layer.The first sandy aquifer of low phase values coincidences with the borehole data.
Figure 3c represents the coverage plot for the model of profile P1 (for location see Fig. 1).According to Günther (2004), the zones of higher coverage values indicate that the model can be reliably derived from the data.Consequently, the maximum depth of sensitive area is about 12 m (Fig. 3c) and beyond this depth the data will lose the capability to resolve the heterogeneity between sedimentary layers.Accordingly, the heterogeneity within the upper aquifer can be well recognized (Fig. 3b).
The inversion results of the other longer 2-D profile P2 showed no great differences to P1.As with the P1 profile, we chose the 0.625 Hz data for presentation and further discussion.The regularization parameter of z weight (w z ) was set to a z /a x = 0.01.The resistivity amplitude |ρ| image exhibits a predominantly layered structure with lateral variations within these layers (Fig. 4a).The phase image of the complex electrical resistivity is shown in Fig. 4b.It exhibits a clear layering in comparison with the resistivity amplitude sections.Based on the coverage plot of the measured data, the maximum depth of the reliably inverted model is 30 m (Fig. 4c).
The high resistivity (> 300 m) complex of 5.6 m thickness (Fig. 4a) can be differentiated into three layers in the phase image (Fig. 4b); a low phase layer (< 4 mrad), which corresponds to sand, between two high phase (> 8 mrad) layers corresponding to peats.Because of the high data coverage at 4.1 m depth (Fig. 4c), the peat layer can be well recognized in comparison with the inversion results of profile P1 (Fig. 3b).The soil layer is followed by a medium resistivity layer of ∼ 18 m thickness, which can correspond to the sandy aquifer.This layer overlays a low resistivity layer (< 40 m) corresponding to Cretaceous marl.The heterogeneity within the sandy aquifers cannot be recognized.On the contrast, the phase shift image (Fig. 4b) represents the medium resistivity layer in the form of two layers; the upper low phase layer (< 4 mrad) corresponding to the first sandy aquifer, which is followed by a high phase layer (> 13 mrad), corresponding to the till layer.The lower sandy aquifer cannot be differentiated from the upper till.This can be attributed to the presence of a high phase layer (till) above the lower aquifer decreasing the resolution with increasing depth and/or to the low data coverage with increasing depth (Fig. 4c).At about 23.6 m depth, this high phase layer is followed by the low phase layer (< 7 mrad) corresponding to the Cretaceous marl.This layer cannot be recognized at the eastern part of the profile, which can be attributed to the bad data coverage (Fig. 4c).

Laboratory measurements and model comparison
In the laboratory, the samples were placed into a transparent sample holder of the SIP Fuchs equipment (Fig. 5).The sample has to be compacted and saturated in the sample  holder, which has to be fixed between two chambers that are filled with water.The pumped water from the drilled borehole was used to carry out the laboratory measurements.The water conductivity was 450 µS cm −1 , which is similar to the groundwater conductivity.Sample holder and chambers were put into a conditioning cabinet in order to guarantee constant temperatures for every measurement (25 increase of ϕ with frequency (f ). Figure 6b shows the IP response (in-phase and quadrature conductivities) at 1.5 Hz of the core samples at different depths.As in Weller et al. (2011) and Weller and Slater (2012), a direct relation between water conductivity and the real component of electrical conductivity was observed.In addition, the imaginary conductivity exhibits a steeper increase at low salinities and asymptotic behavior (flattens) at high water conductivity.Figure 6c and d display the variation of in-phase and imaginary conductivity spectra for the 7.58 m depth core sample at different concentrations of sodium chloride solution.Larger variation of σ is observed with varying salt concentration to more than 115 µS cm −1 (Fig. 6c).A more or less continuous increase of σ with rising frequency is observed (Fig. 6d).Generally, an increase of σ with rising concentration is observed, although the spectral curves of the three highest salt concentrations cannot be clearly separated.As shown in Revil and Florsch (2010) and Revil et al. (2012b), the electrical conductivity data of some selected samples are shown in Fig. 7.The formation factor (F B ) of the samples ranged from 5 to 15 using σ ∞ = 1/F B (σ w + σ ∞ s ) (Revil and Florsch, 2010), where σ ∞ is the high-frequency conductivity and σ ∞ s is the surface conductivity.It can be observed that for core samples of peat layers (Fig. 7, right), where the clay and silt sediments are dominant, the fitting between σ w and the high frequency conductivity σ ∞ is different than other samples (Fig. 7, left).
Figure 8 shows the borehole information (ENG_03) and SIP measurements (field and laboratory data) at the intersection point (x = 154 m, see Fig. 3). Figure 8     till, marl and clay and low phases for the upper sandy aquifer (< 5 mrad).In the 2-D resistivity model, the resistivities of the upper aquifer and the till layer are in the same order of amplitude.In spite of good correspondence of resistivity amplitudes between IP 2-D inversion results and laboratory data, the sounding curve inversion results show slightly lower resistivity values.In both models, thin peat layers at around 2 m depth are detected by high phase values, which coincide with laboratory data.On the other hand, the 2-D phase model shows the peat layer at 4.1 m with higher phase values than the sounding model and laboratory data.Laboratory resistivity data show slightly lower values in the second aquifer than the sounding model and they are slightly higher than 2-D model.For the second aquifer, phase values from laboratory SIP are in general heterogeneous but often slightly higher than the sounding model and lower than the 2-D model.Clearly, a correspondence between the field phase models and laboratory data below till layer, i.e., the second aquifer, cannot be observed.

Hydraulic conductivity (K) estimation
To estimate K from IP results, three different approaches were applied using Eqs.( 5), ( 6) and ( 9).These approaches are meaningful in the saturated zone only (i.e., below 2 m).
For calibration and the sake of clearness, the derivation of K values was focused on the laboratory data for the upper aquifer.Then, the K will be predicted for the second aquifer based on the derived K of the upper aquifer.In order to apply Eq. ( 5), a linear relationship should be achieved with a good correlation coefficient (R 2 ) between the logarithms of aquifer resistivity ρ and K.We calculated the real part of resistivity ρ using Eq. ( 4). Figure 9 shows log-log plot of ρ (at 1.5 Hz) against the measured K using coring sleeves of the nine core samples (green solid circles).While there is a general trend of increasing ρ with increasing K, the relationship is weak (R 2 = 0.4).Consequently, the estimation of K from individual ρ seems to be of limited applicability.
Both Börner and Slater models were applied to predict K from the laboratory data and field inversion results in comparison with the measured Kvalues.Note that the adjustable parameters for each model were derived from the correlation with the measured K values of the first aquifer and then they were applied for the second aquifer.The inversion results of IP field data (Fig. 8) at the intersection point (X = 154 m) were used to calculate K using both models.Equation ( 6) of Börner's model was applied using σ w = 0.015 S m −1 and standard values of a = 0.00475 and l = 0.1, as explained above.The exponent c = 2.8 gave the best least-squares fit for the measured K values below the groundwater level.S por was calculated using the Weller et al. ( 2010b) approach (Eq.8).The second approach for K estimation after Slater and Lesmes (2002) was applied using m = 2.3 × 10 −4 and n = 0.9.
Figure 10 shows the applicability of both Börner and Slater models on the laboratory data in comparison with the measured K values from coring sleeves and grain size analyses.The Börner approach represents bigger variations (∼ 10 −7 to ∼ 10 −2 m s −1 ), while the K estimation after the Slater model is varying in a smaller range (10 −6 to 10 −4 m s −1 ).Focusing on the sandy aquifers, it is clear that great differences can be noticed between the results of both approaches for the first aquifer.For the first aquifer, the calculated K values from Börner's model are, in general, higher than the measured K values.On the contrast, a good correspondence   both sounding and profiling data are likely to be caused by the Cretaceous marl.
Based on the above results, it is clear that the use of values from a CPA model, i.e., single frequency data as representative of the spectral response, is not valid to predict K. Consequently, the spectra were investigated using the DD to calculate the following parameters: DC resistivity ρ , total chargeability (M) and different representative relaxation times (i.e., maximum relaxation time τ max ; median relaxation time τ 50 and logarithmically weighted average relaxation time τ lw τ lw =  shape range of complex conductivity spectra and was stable for all samples and all salinities.Figure 12a shows an example of raw data (amplitude and phase) from a core sample at 6.68 m depth.Figure 12b shows the correction of phase values after electromagnetic removal.In addition, Fig. 12d reflects a variation in spectral chargeability with τ by using DD fitting (Fig. 12c).For this sample, the minimum and maximum values of τ are 10 −4 s and 1 s, respectively, and the smoothness factor (λ) used is 150.Phase data were quite accurate with stacking errors of about 0.1 mrad, which allowed the spectral analysis of very small residual phases (Fig. 12).Due to the number of frequencies measured, the derived relaxation time spectrum is relatively independent of statistical errors.
For our samples of the first aquifer, the relation of τ (i.e., τ max , τ 50 and τ lw ) and the measured K was studied.While there is a general trend of increasing K with τ , the relationships are weak.Figure 13 shows an example of the τ -K relationship using different values of τ .Although the τ max -K and τ lw -K relations show a better correlation than τ 50 -K, the correlation coefficients are still low (R 2 = 0.3 and 0.22, respectively).Correspondingly, the use of individual τ to predict K using a power law relation will not be applicable.Moreover, the total of nine measured K values is not enough to apply the multivariate power law relation after Weller et al. (2010a).

A new approach
In general, the relaxation time τ is the time required to establish a stationary ionic distribution in the Stern layer coating a grain of diameter d under the action of a static electrical field (Revil and Florsch, 2010).The relation between the relaxation time with the square of the grain size is represented in Eq. ( 1), which was confirmed by Titov et al. (2002) and Leroy et al. (2008).
The hydraulic and electrical conductivity of an aquifer depend on several factors; such as pore-size distribution, grain size distribution and fluid salinity.If we consider the ρ as a measure of grain size and pore-size distribution (e.g., Kresnic, 2007;Chandra et al., 2008;Attwa et al., 2009;and Attwa, 2012), the aquifer resistivity could be well correlated with the K values of the aquifer (Eq.5).Similarly, if we consider τ as a measure of the pore radius and pore size characteristics (Eqs. 1, 2), a link between τ and K can be considered.The Debye decomposition of IP spectra yields a characteristic relaxation time τ .These can be expressed in the form of Kα(ρ ) B (τ ) C .Accordingly, the K can be calculated using where the A coefficient and the B and C exponents are determined by a regression analysis of the logarithmic quantities τ and K.This methodology was adopted for the estimation of the K when there is a general trend of increasing or decreasing both individual ρ and τ with K.

Correlation studies and testing methodology
In the cases of ρ -K or τ -K relationships, the correlations are weak but the there is a general increase of ρ and τ with the measured K values (Figs. 9 and 13, respectively).Equation (13) was applied to the lab data collected from the upper (first) aquifer.Then, the coefficients A and the exponents (B and C) were used to predict the K values of the lower aquifer (i.e., the second one).The ρ τ -K relation was studied using τ max , τ 50 (corresponds to the mean relaxation time as proposed by Zisser et al. (2010b) and τ lw .A strong correlation coefficient (R 2 = 0.81) was achieved between ρ τ lw and the measured K (Fig. 14) for nine samples of the first aquifer.A multivariate power law relationship (A(ρ ) B (τ lw ) C ) was

Fig. 14.
Log-log plot of τ max , τ 50 and τ lw multiplied by ρ against the measured K of nine core samples from the first aquifer (see Fig. 1).
examined and it was observed that the exponents B and C are nearly equal; K = 0.00012 (ρ ) 1.08 (τ lw ) 0.8 .However, due to the limited number of data points we fixed them to be equal (i.e., K = A(ρ τ ) B ), thus decreasing the number of unknowns.A and B were calculated for the upper aquifer and they were 0.0006 and 0.73, respectively.Based on these results, Eq. ( 13) was applied for the all core samples of both aquifers using the deduced coefficient A and exponent B.
Figure 10 shows the calculated K values using Eq. ( 13) in comparison with the measured K values and the calculated K values from single frequency approaches.In comparison with Eqs. ( 6) and ( 9), Eq. ( 13) shows the lowest mean-squared logarithmic error of fit (Tong and Tao, 2008), where N is the number of samples (in this paper it is 29), K m (i) and K c (i) are the measured and computed K by conventional measurements and various approaches, respectively, and δ is the error factor.For Börner and Slater models, the δ values were 3.5 and 3.6, respectively.On the other hand, the δ was 1.4 for the proposed approach (Eq.13).Accordingly, Eq. ( 13) is at least applicable at this test site (Fig. 10).

Discussion
Since it is hard to gain a full control of the subsurface structures in natural geological environments during the 2-D inversion, the discrepancy between 1-D and 2-D inversion results can be observed (Fig. 8).Because the 2-D inversion process is carried out using a fixed mesh size, the small IP response can be embedded within the general IP response of the layer, which could be amplified during the inversion process.On the other hand, the 1-D inversion shows the IP response for a part of a layer with depth and accordingly, the lateral effect of IP response will be ignored.
Based on the inversion results of both sounding and profile P1, it is clear that the upper aquifer can be well defined, however we still have limitations in imaging the lower one.These limitations include an electromagnetic coupling with increasing current electrode spacing of the sounding and the penetration depth of the 2-D imaging.In general, our field inversion results indicate that the upper and lower boundaries of the upper (first) aquifer can be well defined from the phase values.The upper boundary of the lower sandy aquifer (i.e., the second one) appears deeper than in the borehole, which can be attributed to the observed misfit between the raw and fitted phase angle data (Fig. 2, right).Also, there is a difference of the phase behavior of the peat layer between the laboratory, sounding and profiling results (Fig. 8, right).Zanetti et al. (2011) and Ponziani et al. (2011) reported high phase values for buried tree root samples using laboratory experiments.They showed that a regular distribution of the pore size in diffuse porous woods leads to a stronger polarization effect.
As in Weller et al. (2006) and Revil et al. (2012b), a considerable variation in the spectra of conductivity amplitude and phase shift can be observed in laboratory investigations of peat, wood and soil samples.Weller et al. (2006) concluded that characteristic features of phase spectra can be used to detect wooden relics of trackways or pile dwelling from laboratory investigations.Revil et al. (2012b) developed a quantitative model to investigate the frequency domain induced polarization response of suspensions of bacteria and bacteria growth in porous media.They showed that the induced polarization of bacteria (α polarization) is related to the properties of the electrical double layer of the bacteria.The high IP effect can be attributed to the surface of the bacteria, which is highly charged due to the presence of various structures extending away from the membrane into the pore water solution.Consequently, the dead cells (i.e., organic material) will show low IP effects.
The inspection of both sensitivity models and the inversion results of IP data is efficient in evaluating the reliability of the interpretation.The IP profiles (Figs. 3,4) show a high polarization effect of the peat layer compared to laboratory and sounding results.The coverage model demonstrates that the reliable results concentrate up to a maximum depth of 16.5 m, which is only sufficient to cover the upper sandy aquifer.Because we have bad data coverage and electromagnetic coupling effects with increasing electrode spacing, the second aquifer cannot be imaged well.Accordingly, accurate K values of the lower aquifer from the IP field data cannot be expected.Figure 8 shows slight differences between field data inversion results and laboratory data up to the till layer (∼ 13 m depth).Koch et al. (2011) stated that changes in compaction and sorting of samples cause a certain shift in the SIP response but did not fundamentally alter the overall picture.Accordingly, the modification of pore space and grain size characteristics of original samples through compaction of preparing the lab sample will change the SIP response, while the chemistry of the saturation pore fluid is kept constant.
Regarding the above results, the proposed IP-K relationships based on single frequency models appear inappropriate for the unconsolidated heterogeneous sediments under investigation.Clearly, there are limitations of using both Börner and Slater models.In order to apply the single frequency models (Eqs.6, 9) of both Börner et al. (1996) and Slater and Lesmes (2002), certain conditions are required: (a) a constant phase angle over a wide frequency range, (b) a direct relation between S por and σ and (c) an inverse relation between the measured K and σ .Figure 6 shows noticeable differences in phase angle with frequencies and it would appear inappropriate to use a single frequency complex measure as a representative of IP spectra.The correlation between S por and σ cannot be proved because S por values were not measured (i.e., by BET measurements).Here, the correlation between S por and σ cannot be proved because S por values were not measured (i.e., by BET measurements).Because Eq. ( 8) after Weller et al. (2010b) was derived from an extensive sample database, the direct relation can be shown between S por and σ for unconsolidated sediments.Moreover, Fig. 9 shows that the relationship between σ and K is weak and not linear, although a general decrease of σ with K can be observed.The existence of laboratory data over one decade could be the reason for this weak relation.Therefore, the models of Börner and Slater seem to be of limited applicability at this test site.Compared with the overall lithology and measured K values, the results from the Slater model appear more realistic than those from the Börner model, but, in general, the predicted K values are too small.Since phase spectra of unconsolidated sediments typically show frequency dependence, it is difficult to obtain a constant phase angle behavior over a limited frequency range for inhomogeneous unconsolidated aquifers; multiple frequency data are required to study the spectral behavior and to derive the K values.Our laboratory measurements of inhomogeneous sandy aquifers show a weak relationship between the measured hydraulic conductivity and the maximum, median and weighted average relaxation time.Independently of the soil texture, soil compaction has an influence on both the conductive properties (Seladji et al., 2010) and the polarization response (Koch et al., 2011) of porous media.Accordingly, the weak correlation between K and τ can be attributed to the compaction and sieving of the core samples from preparing the laboratory measurements (Koch et al., 2011).Consequently, other physical parameters, which show a relation with the measured K, are required to support the K-τ relationship.Since the aquifer resistivities show a positive direct relation with the measured K values, a multiplication of τ and aquifer resistivities shows a good power law relationship to estimate the K.In our study, the logarithmically weighted average relaxation time shows superior results of predicting the K compared with the mean relaxation time used by Weller et al. (2010a) and the median relaxation time recommended by Zisser et al. (2010b).

Conclusions
The goal of this study was to use SIP field and laboratory measurements to characterize the hydrogeological conditions of the Schillerslage test site in Germany.The SIP method provided a tool for imaging subsurface hydrogeological parameters of unconsolidated inhomogeneous sediments.The improved Debye decomposition procedure was used to calculate (τ ) values from the measured spectra.We have reached the following conclusions.
1.In our particular case, the IP inversion results indicated that the understanding of SIP signatures from laboratory data can be achieved in field applications.The results here suggested that SIP is clearly able to differentiate lithological heterogeneities within unconsolidated sediments.
2. The combination of field and laboratory IP measurements was recommended to evaluate the geological situation and the accuracy of the interpretation (e.g., changes in compaction and sorting of samples resulted in a certain shift in the SIP response between laboratory and field inversion results).The mentioned recommendation provided a means to (i) overcome the spectral limitations resulting from EM coupling effects, which mask the SIP response at higher frequencies, and (ii) to examine the influence of physical and chemical properties of laboratory samples on SIP response.
3. Hydraulic conductivity estimates showed that single frequency models (Börner et al., 1996, andSlater andLesmes, 2002) used to predict K were inappropriate according to our investigations.However, the Slater and Lesmes (2002) model yielded more realistic results than the model of Börner et al. (1996).These can be attributed to (i) the different type of sediments used to derive the different approaches and (ii) because the assumption that the phase shift is constant over a wide frequency range was not achieved in our data.Therefore, the prediction of K was examined using the relaxation time (τ ).It was observed that a direct relation between the K and τ was not successful.In addition, the relation between the real resistivity part (ρ ) and K was weak.Another practical application of ρ τ -K was tested.
4. Our results revealed that the integration of the logarithmic weighted average relaxation time (τ lw ) and DC resistivity yielded a power law relation to predict the K of unconsolidated sandy aquifers particularly if only few data points are present for correlation.However, it would be more meaningful if the above relation is tested in areas with diverse geological environments.Moreover, with the number of samples increasing, a multivariate power law relation (K = A(ρ ) B (τ lw ) C ) is expected to improve the prediction of K values and, consequently, this relation should be examined.
5. The present approach offers the possibility to use the field data to monitor hydrogeological parameters using the SIP components.A limited range of hydraulic conductivity variation over the study area can be observed and, consequently, the derived equation in this study is not expected to apply to other areas but the methodology will.
6. Further work including lab measurements is required to study the spectral behavior of the peat layer.In addition, application of recent universal models for calculating the specific surface area and the permeability using the formation factors are desirable in comparison with the approaches we used in this paper.Taking these points into account can improve the accuracy of predicting the hydraulic conductivity measurements.
In this framework a further use of the recent models (universal functions) can be suggested to account for the dependence of the in-phase and quadrature conductivities on the porosity, cation exchange capacity (or specific surface area of the material), and salinity of the pore water for simple supporting electrolytes like NaCl in isothermal conditions (Revil, 2013;Revil et al., 2013a, b).

Fig. 2 .
Fig. 2. SIP sounding inversion results.Apparent resistivity amplitude (a) and phase angle (b) data and the inversion results (c and d, respectively) in comparison with well ENG_03 (see Fig. 1).

Fig. 5 .
Fig. 5. Layout of laboratory SIP measurements showing how the sample holder is measured using the SIP Fuchs equipment.

Fig. 5 .
Fig. 5. Layout of laboratory SIP measurements showing how the sample holder is measured using the SIP Fuchs equipment.

Fig. 6 .
Fig. 6.(a) Exemplary SIP response at three samples at different depths (2.8 m, 8.82 m and 11.8 m) and (b) in-phase and quadrature conductivities (at 1.5 Hz) at different depths.(c and d) Inphase and quadrature conductivities at different water salinities of the collected sample at 7.58 m depth.

Fig. 6 .
Fig. 6.(a) Exemplary SIP response at three samples at different depths (2.8, 8.82 and 11.8 m) and (b) in-phase and quadrature conductivities (at 1.5 Hz) at different depths.(c, d) In-phase and quadrature conductivities at different water salinities of the collected sample at 7.58 m depth.

Fig. 10 .
Fig. 10.Hydraulic conductivity of samples from borehole Eng_03.Black line: own approach.Green line: from coring sleeves and grain size analysis using Kozeny equation.Pink line: from the Slater and Lesmes (2002) model.Blue line: from the Börner et al. (1996) model.

Fig. 12 .
Fig. 12. SIP lab data of one core sample at 6.68 m depth.(a) Resistivity amplitude (|ρ|).(b) Measured phase (φ) (blue line) values, electromagnetic (EM) term (green line) and phase values (red line) after removing EM effect.(c) The fitting between the phase values after removing EM effect (blue line) and the fitted curve (red line).(d) The variation of spectral chargeability with relaxation time by using DD model.

Fig. 12 .
Fig. 12. SIP lab data of one core sample at 6.68 m depth.(a) Resistivity amplitude (|ρ|).(b) Measured phase (ϕ) (blue line) values, electromagnetic (EM) term (green line) and phase values (red line) after removing EM effect.(c) The fitting between the phase values after removing EM effect (blue line) and the fitted curve (red line).(d) The variation of spectral chargeability with relaxation time by using DD model.

Fig. 13 .Fig. 14 .
Fig.13.Log-log plot of the maximum relaxation time (τ max ), median relaxation time (τ 50 ) and weighted average relaxation time (τ lw ) against the measured K of nine core samples from the first aquifer (see Fig.1).