# Subsurface imaging of water electrical conductivity, hydraulic permeability and lithology at contaminated sites by induced polarization

Subsurface imaging of water electrical conductivity, hydraulic permeability and lithology at... Summary At contaminated sites, knowledge about geology and hydraulic properties of the subsurface and extent of the contamination is needed for assessing the risk and for designing potential site remediation. In this study, we have developed a new approach for characterizing contaminated sites through time-domain spectral induced polarization. The new approach is based on: (1) spectral inversion of the induced polarization data through a reparametrization of the Cole–Cole model, which disentangles the electrolytic bulk conductivity from the surface conductivity for delineating the contamination plume; (2) estimation of hydraulic permeability directly from the inverted parameters using a laboratory-derived empirical equation without any calibration; (3) the use of the geophysical imaging results for supporting the geological modelling and planning of drilling campaigns. The new approach was tested on a data set from the Grindsted stream (Denmark), where contaminated groundwater from a factory site discharges to the stream. Two overlapping areas were covered with seven parallel 2-D profiles each, one large area of 410 m × 90 m (5 m electrode spacing) and one detailed area of 126 m × 42 m (2 m electrode spacing). The geophysical results were complemented and validated by an extensive set of hydrologic and geologic information, including 94 estimates of hydraulic permeability obtained from slug tests and grain size analyses, 89 measurements of water electrical conductivity in groundwater, and four geological logs. On average the IP-derived and measured permeability values agreed within one order of magnitude, except for those close to boundaries between lithological layers (e.g. between sand and clay), where mismatches occurred due to the lack of vertical resolution in the geophysical imaging. An average formation factor was estimated from the correlation between the imaged bulk conductivity values and the water conductivity values measured in groundwater, in order to convert the imaging results from bulk conductivity to water conductivity. The geophysical models were actively used for supporting the geological modelling and the imaging of hydraulic permeability and water conductivity allowed for a better discrimination of the clay/lignite lithology from the pore water conductivity. Furthermore, high water electrical conductivity values were found in a deep confined aquifer, which is separated by a low-permeability clay layer from a shallow aquifer. No contamination was expected in this part of the confined aquifer, and confirmation wells were drilled in the zone of increased water electrical conductivity derived from the geophysical results. Water samples from the new wells showed elevated concentrations of inorganic compounds responsible for the increased water electrical conductivity in the confined aquifer and high concentrations of xenobiotic organic contaminants such as chlorinated ethenes, sulfonamides and barbiturates. Electrical properties, Permeability and porosity, Hydrogeophysics, Inverse theory, Tomography 1 INTRODUCTION Contamination of groundwater and surface water by heavy metals, nutrients or xenobiotic organic compounds in urban areas is a common problem all over the world (Panagos et al.2013). In particular, contaminated sites (e.g. former industrial facilities or dump sites, old landfills, gasoline stations or dry cleaning facilities) can generate groundwater plumes posing a risk to water resources (Basu et al.2006; Ellis & Rivett 2007; Bjerg et al.2011). For site investigations, risk assessment and adequate remedial design, a characterization of the contaminated sites and plumes is needed (Sims 1990; Barcelona 1994). This requires information on the spatial distribution of the contaminant plume, local geology and hydrogeological properties (e.g. permeability, k), which are conventionally estimated from analysis of groundwater samples, of soil samples and of aquifer tests (EPA 1991; Cameron 1992; Benson & Yuhr 2016a). These approaches require many drillings to sufficiently characterize field sites, which is often unfeasible at large sites (Benson & Yuhr 2016b). As an alternative approach, non-invasive geophysical methods can be used in combination with drillings for improved 3-D characterization of contaminant plumes. Among the various geophysical techniques, the direct current (DC) resistivity method has been used for mapping groundwater contamination. The presence of contamination affects the formation resistivity depending on the type and concentration of contaminants. The increase in resistivity is usually associated with the presence of mobile (pools) or residual non-aqueous phase liquids (NAPLs; Deryck et al.1993; Yang et al.2007; Johansson et al.2015). A decrease in the resistivity is typically associated with ionic compounds, which may be linked to contamination and/or to different redox or biodegradation processes (Atekwana et al.2005; Chambers et al.2006; Junejo et al.2015; Maurya et al.2017). The DC resistivity method suffers from well-known limitations, that is, its inability to discriminate between low-resistivity formations (e.g. clay) and high salinity of groundwater. With the IP method, which measures the capacitive nature of the subsurface (Binley 2015), additional information about the surface conductivity is gained (to avoid confusion, here and throughout the manuscript, ‘conductivity’ always refers to electrical conductivity, while ‘permeability’ is used when referring to hydraulic properties). This helps to discriminate lithology-driven resistivity variations from variations driven by pore water. In particular, Weller et al. (2013) identified a strong linear relationship between the real part of surface conductivity and the imaginary conductivity, and discovered how this relationship can be used to improve the estimation of the true formation factor and the groundwater conductivity when an IP measurement is also made. Over the past two decades, the spectral nature of the IP response has increasingly been used as an exploration tool in environmental and hydrogeological investigations (see Kemna et al.2012, for an overview). Spectral IP (SIP) signals can be measured in both time domain (TD) and frequency domain (FD). In TD, SIP measurements are usually combined with DC resistivity measurements using similar field procedures and instrumentation as when measuring DC data alone. Several studies about inversion and processing of TD SIP data have recently been published (Hördt et al.2006; Hönig & Tezkan 2007; Fiandaca et al.2012; Fiandaca et al.2013; Olsson et al.2016; Fiandaca et al.2017a) and successful case studies have been presented where the TD SIP method was applied for mapping lithology (Gazoty et al.2012b; Johansson et al.2016), landfill-waste materials (Gazoty et al.2012a) and contaminated sites (Johansson et al.2015; Sparrenbom et al.2017). In addition to complementing the DC method for lithological discrimination and characterization of contaminated sites, the IP method has been used to estimate hydraulic permeability in the laboratory (e.g. Binley et al.2005; Slater 2007; Revil & Florsch 2010; Zisser et al.2010; Weller et al.2015) and in the field (e.g. Kemna et al.2004; Hördt et al.2007; Hördt et al.2009; Attwa & Günther 2013; Binley et al.2016). In particular, Weller et al. (2015) proposed empirical equations for permeability estimation derived from an extensive set of samples for consolidated and unconsolidated sediments. The applicability of the proposed equation for unconsolidated sediments was verified in the field, in the presence of heterogeneity in both lithology and water chemistry, at the same site presented in this study with the El-log technique (Fiandaca et al.2017b), that is, a borehole technique for acquiring logging-while-drilling DC resistivity, TD spectral IP (SIP) and gamma radiation data (Sørensen & Larsen 1999; Gazoty et al.2012a). The high vertical resolution of the El-log technique matched the lithological variability at the site, minimizing the ambiguity in the result interpretation due to lack of geophysical resolution. A very good correlation within on average one order of magnitude was found between the IP-derived permeability estimates and those derived using grain size analyses and slug-tests, with similar depth-trends and permeability contrasts. These new results on the use of TD SIP data for characterization of hydraulic permeability, the relation between real and imaginary surface conductivity derived in the laboratory, and the recent improvements in processing and inversion of TD SIP data paved the way for the development of a new approach for the characterization of contaminated sites through TD SIP. The new approach is based on: (1) spectral inversion of the TD SIP data through a new reparametrization of the Cole–Cole model (Cole & Cole 1941; Pelton et al.1978), which disentangles the electrolytic bulk conductivity from the surface conductivity for delineating the (inorganic) contamination plume; (2) using the relationship derived from Weller et al. (2015) for estimating the hydraulic permeability in the field directly from the inversion parameters, without any calibration; (3) the use of geophysical imaging results for supporting the geological modelling. The new approach was tested on a TD SIP survey close to the Grindsted stream (Denmark), where the contaminated groundwater from a factory site discharges to the stream (Sonne et al.2017; Rønde et al.2017). The new approach was complemented and validated by comparing the geophysical results with an extensive set of hydrologic and geologic information. In particular, the geophysical results were compared with 94 estimates of hydraulic permeability obtained from slug tests and grain size analyses, 89 measurements of water conductivity in groundwater and four geological logs. 2 METHODOLOGY 2.1 Induced polarization and hydraulic permeability In the induced polarization method the low-frequency capacitive properties of the earth are measured as a frequency-dependent complex resistivity (in FD) or a decay response (in TD) when the medium is excited by a time-varying electric current. In a porous medium free of metallic particles, the electrical conduction takes place through the fluid that fills the interconnected pore spaces (electrolytic bulk conduction) and through the electrical double layer (surface conduction), and these conductivities are usually assumed to add in parallel into the total complex conductivity σ* (e.g. Lesmes & Frye 2001):   \begin{eqnarray} {\sigma ^*} ( f ) &=& {\sigma _{{\rm{bulk}}}}( {{\sigma _{w}}}) + \sigma _{{\rm{surf}}}^* ( {{\sigma _{w},f}}) \nonumber\\ &=& \frac{{{\sigma_{w}}}}{F}\ + \sigma _{{\rm{surf}}}^{'}\left( {{\sigma _{w},f}} \right) + i\sigma _{{\rm{surf}}}^{''}\left( {{\sigma _{w},f}} \right) \end{eqnarray} (1)where the bulk conduction σbulk is expressed by the water conductivity σw divided by the formation factor F and the dependence of the complex surface conductivity $$\sigma _{\rm {surf}}^*$$ over the frequency f and the water conductivity σw is stated expressly. The real and imaginary components of the surface conductivity are not independent: using a database composed of 63 sandstone and unconsolidated sediment samples covering nine independent investigations, Weller et al. (2013) identified a strong linear relationship between $$\sigma _{\rm {surf}}^{'}$$ determined from multisalinity resistivity measurements and $$\sigma _{\rm {surf}}^{''}$$ measured with IP at a frequency of about 1 Hz:   $$\ {\sigma ''_{\rm {surf}}} = \ l \cdot {\sigma '_{\rm {surf}}}$$ (2)with l = 0.042 ± 0.022 (dimension less). Furthermore, Weller et al. (2013) found that the salinity dependency of the real surface conductivity parallels the salinity dependency of the imaginary surface conductivity, which can be expressed as (Weller et al.2011; Weller et al.2015):   $${\sigma ''_{{\rm{surf}}}} ( {{\sigma _w}}) = {\sigma ''_{{\rm{surf}}}}( {{\sigma _f}}) \cdot \frac{1}{{{C_f}}}\sqrt {\frac{{{\sigma _w}}}{{{\sigma _f}}}}$$ (3)where σf represents the salinity of a reference NaCl solution and Cf accounts for the possible differences in ionic species in the reference and actual solution. IP-based permeability prediction methods are based either on the imaginary conductivity $$\sigma _{\rm {surf}}^{''}$$ or on relaxation time τ (a quantity used to represent the characteristic hydraulic length scale e.g. Revil (2012). The premise of estimating k using $$\sigma _{\rm {surf}}^{''}$$ is based on its strong relationship with surface area normalized to the pore volume (Spor), which holds the fundamental basis for derived empirical relationships between k and induced polarization in laboratory studies (Börner et al.1996; Slater & Lesmes 2002a; Weller et al.2015). Recently, Weller et al. (2015) investigated a database consisting of 114 globally collected samples. They avoided the indirect k-estimation through Spor and suggested direct correlations between k and the electrical parameters. For unconsolidated sediments, they proposed the following empirical equation:   $$k\ = \frac{{1.08 \times {{10}^{ - 13}}}}{{{F^{1.12}} \cdot {{\left( {{{\sigma ''}_{{\rm{surf}}}}\left( {{\sigma _f}} \right)} \right)}^{2.27\ }}\ }}.$$ (4)Substituting eq. (3) and using the Archie's law $$({{\sigma _{{\rm{bulk}}}} = \frac{{{\sigma _W}}}{F}})$$, we can re-write eq. (4) as:   $$k\ = \frac{{5.80 \cdot {{10}^{ - 16}}}}{{{C_f}^{2.27}}} \cdot \frac{{{\sigma _{{\rm{bulk}}}}^{1.12}}}{{\sigma {{_{{\rm{surf}}}^{''}}^{2.27\ }}}} \cdot {\sigma _w}^{0.015}$$ (5)where in the last equality the explicit dependence of $$\sigma _{\rm {surf}}^{''}$$ on σw was omitted and the value of σf = 100 mS m − 1 was used. The permeability estimation through eq. (5) depends weakly on water conductivity: a 10-fold and 100-fold variations in water conductivity cause only approximately 3.5 per cent and 7 per cent variations in permeability, respectively. For this reason, in the following the spatial variability of σw is disregarded in permeability computation, and the uniform value σw = 47 mS m − 1 is used, that is, the average value at the site (the range at the site is 11–86 mS m − 1). Weller et al. (2011) suggest Cf = 2 for CaCl2 and Cf = 1 for NaCl; for other ions, no suggestion was made. Considering that a mixture of cations and anions are present in the field-collected water samples with varying molecular concentration, it is difficult to apply an appropriate correction. Therefore, in our k-estimation the value Cf = 1 is used regardless of the chemical composition of the water. 2.2 Parametrization of induced polarization The spectral variation of the complex conductivity is often parametrized, and the Cole–Cole model (Cole & Cole 1941; Pelton et al.1978) is often used for modelling IP data acquired in the field (Fiandaca et al.2013; Kemna et al.2014; Günther & Martin 2016). The Cole–Cole model in its conductivity form is defined as (e.g. Tarasov & Titov 2013)   $${\sigma ^{\rm{*}}} (f) = {\sigma _0}\ \cdot \left[ {1 + \frac{{{m_0}}}{{1 - {m_0}}} \cdot \left( {1 - \frac{1}{{1 + {{\left( {i2\pi f{\tau _\sigma }} \right)}^C}}}} \right)} \right]$$ (6)where σ* is the complex conductivity, σ0 is the DC conductivity, m0 is the intrinsic chargeability, τσ is the relaxation time, C is the frequency exponent, f is the frequency and i is the imaginary unit. Using Monte Carlo Markov Chain analysis, Fiandaca et al. (2017a) have shown that the parameters of the Cole–Cole model are strongly correlated when inverting IP data (in particular m0 and C) and that smaller parameter correlations and better resolution can be achieved (in both FD and TD) when re-parametrizing the Cole–Cole model, replacing the parameter m0 with the maximum imaginary conductivity (MIC) $$\sigma _{{\rm{max}}}^{{\rm{''}}}$$ of the Cole–Cole spectrum (Fig. 1b). Figure 1. View largeDownload slide Spectrum of the BIC model computed using $${\sigma _{{\rm{bulk}}}} = {\rm{\ }}2{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ \ }}\sigma _{{\rm{max}}}^{{\rm{''}}} = {\rm{\ }}0.5{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ }}{\tau _\sigma } = 0.05{\rm{\ s\ and\ }}C = 0.5$$. (a) Real conductivity σ'(black curve), obtained as the sum of the bulk conductivity σbulk (magenta curve) and the surface real conductivity$${\rm{\ }}\sigma _{{\rm{surf}}}^{\rm{'}}$$. The water conductivity value σw with formation factor F = 5 is also shown. (b) Imaginary conductivity σ''. Figure 1. View largeDownload slide Spectrum of the BIC model computed using $${\sigma _{{\rm{bulk}}}} = {\rm{\ }}2{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ \ }}\sigma _{{\rm{max}}}^{{\rm{''}}} = {\rm{\ }}0.5{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ }}{\tau _\sigma } = 0.05{\rm{\ s\ and\ }}C = 0.5$$. (a) Real conductivity σ'(black curve), obtained as the sum of the bulk conductivity σbulk (magenta curve) and the surface real conductivity$${\rm{\ }}\sigma _{{\rm{surf}}}^{\rm{'}}$$. The water conductivity value σw with formation factor F = 5 is also shown. (b) Imaginary conductivity σ''. The resulting model, named the MIC model, is defined in terms of the parameters mMIC:   $${{\boldsymbol{m}}_{{{\bf MIC}}}} = \{ {{\sigma _0},\sigma _{{\rm{max}}}^{''},\ {\tau _\sigma },\ C} \}$$ (7) In eq. (7) the parameter σ0 represents the total DC conductivity, that is, the sum of the bulk conductivity σbulk and the DC surface conductivity $$\sigma _{{\rm{surf}}}^{'}( {f\ = \ 0} )$$ (see eq. 1), so when imaging the DC conductivity (or, equivalently, DC resistivity), it is difficult to discriminate between low-resistive formations (e.g. clay) and increased salinity of the groundwater. For this reason in this study, we introduce a new Cole–Cole reparametrization, that is, the bulk and (maximum) imaginary conductivity (BIC) model, defined in terms of the parameters mBIC  $${{\boldsymbol{m}}_{{{\bf BIC}}}} = \{ {{\sigma _{{\rm{bulk}}}},\sigma _{{\rm{max}}}^{''},\ {\tau _\sigma },\ C} \}$$ (8)where the assumption is that real and imaginary surface conductivity of eq. (1) are proportional through the proportionality factor l expressed in eq. (2) and that the frequency dependence of the complex conductivity obeys the Cole–Cole model. In Weller et al. (2013), the proportionality factor is found between the DC surface conductivity and the imaginary conductivity at 1 Hz, but the variability of the proportionality with frequency is not discussed. Considering that the imaginary conductivity of the Cole–Cole model reaches a maximum at the frequency f = 1/2πτσ, we decided to enforce in the BIC model the proportionality between the real and imaginary surface conductivity at this frequency. This is a conservative choice: in this way the ratio between the surface imaginary conductivity and real conductivity never exceeds the factor l of eq. (2). On the contrary, enforcing the proportionality at f = 1 Hz would imply a ratio well above l at the peak frequency f = 1/2πτσ for models with τσ ≫ 1. For any given set of BIC model, the corresponding conductivity Cole–Cole parameters can be easy retrieved from the BIC parameters. First, let's define a and b as   $$a\ = \ - {\rm{imag}}\left( {\frac{1}{{1 + {i^C}}}} \right)$$ (9)  $$b\ = \frac{{{m_0}}}{{1 - {m_0}}}\$$ (10) Considering that $${\rm{real}}\ ( {\frac{1}{{1 + {i^C}}}} ) = \ 0.5$$ regardless of the C value, eq. (6) can be written at f = 1/2πτσ as   $${\sigma ^*}( {f\ = 1/2\pi {\tau _\sigma }}) = {\sigma _0}\ \cdot \left[ {1 + b \cdot \left( {0.5 + i \cdot a} \right)} \right]$$ (11) Enforcing the proportionality of eq. (2) at f = 1/2πτσ  and considering that by definition $$\sigma _{{\rm{max}}}^{''} = \ {\rm{imag}}( {{\sigma ^*}( {f\ = 1/2\pi {\tau _\sigma }} )} )$$, we obtain from eq. (1):   $${\sigma _0} \cdot \left[ {1 + b \cdot \left( {0.5 + i \cdot a} \right)} \right]\ = {\sigma _{{\rm{bulk}}}}\ + \sigma _{{\rm{max}}}^{''}/l + \ i \cdot \sigma _{{\rm{max}}}^{''}.$$ (12) The imaginary part of eq. (12) can be solved for σ0:   $${\sigma _0} = \frac{\sigma _{\rm{max}}^{''}} {a \cdot b}.$$ (13) Substituting eq. (13) into eq. (12), it is possible to solve for b from the real part of the resulting equation:   $$b = \frac{l \cdot \sigma _{\rm{max}}^{''}/a} {{\sigma _{\rm{max}}^{''}} + l \cdot {\sigma _{\rm{bulk}}} - 0.5 \cdot l \cdot \sigma _{\rm{max}}^{''}/a}$$ (14) Finally, m0 can be retrieved as a function of b from eq. (10):   $${m_0} = \frac{b}{{1 + b}}\$$ (15) For example, Fig. 1 shows the spectrum of the BIC model defined by the parameter values $$\{ \sigma _{\rm{bulk}} = 2\ {\rm{mS}}\ {\rm{m}}^{ - 1},\ \sigma _{\rm{max}}^{''} = 0.5\ {\rm{mS}}\ {\rm{m}}^{ - 1},\ \tau _\sigma = 0.05\ {\rm{s}}, C = 0.5 \}$$. The σ0 and m0 values of the corresponding Cole–Cole model are σ0 = 12.7 mS m − 1 and m0 = 160 mV V − 1, and the DC conductivity is almost completely dominated by the surface conductivity. 2.3 Time-domain spectral induced polarization TD spectral IP is an extension of the DC resistivity method in which the capacity of the ground is measured using square current pulses. A TD SIP signal can either be measured as rising of the potential during the current on time or decaying of potential after the current is turned off. The measurement using the former approach is referred as 100 per cent duty cycle measurement and the latter as 50 per cent duty cycle measurement (Olsson et al.2015). The possibility of extracting the spectral information contained in the TD signal depends considerably on the number of decades acquired (Madsen et al.2016, 2018), and nowadays commercial instruments exist that sample the full-waveform potential and current signals at high sampling rate (typically in the kHz range) for the entire measurement time, allowing for advanced signal processing. The recorded full-waveform data often contain harmonic noise from the power distribution grid, spikes from, for example, animal fences and self-potential background drift. Harmonic noise makes it impossible to measure earlier than 20 ms (in a 50 Hz environment) and self-potential drift distorts the signal at later times: this limits the extraction of spectral information from the signal. In this study full-waveform data sampled at 3750 Hz were acquired and processed following Olsson et al. (2016) for removal of harmonic noise, spikes, self-potential drift and tapered windowing, in order to retrieve unbiased TD SIP data from a few milliseconds after current switch. Apparent DC resistivity values and full-decay gated curves (i.e. apparent chargeability values with time) represent the data space for the inversion of TD SIP data. 2.4 2-D inversion and imaging of hydraulic permeability The inversion procedure is carried out using the code by Auken et al. (2015), where the apparent resistivity values and the IP full-decays are inverted simultaneously in 2-D following Fiandaca et al. (2013), with an accurate description of transmitter waveform and the receiver transfer function for an unbiased estimation of the spectral parameters (Fiandaca et al.2012). Along the 2-D section, the model space is defined cell by cell in terms of a parametrization of the FD complex conductivity (The MIC and BIC parametrizations are the ones used in this study). The same vertical model discretization (with an increased number of layers in the big-scale profiles), model constraints and starting model were used for all the inversions, as well as the same stopping criterion for the iterative inversion (i.e. a relative variation of the objective function between consecutive iterations below 2 per cent). Furthermore, the same model for the data standard deviations (STD) has been adopted, that is, 1 per cent on apparent resistivity and 10 per cent on chargeability, plus a noise floor of 0.1 mV (Olsson et al.2015). L2 smoothness vertical/horizontal constraints were applied, with values 2.0 and 1.2, respectively. The constraint values represent the relative vertical/lateral variation of the parameters that weights the roughness misfit in the objective function (Auken et al.2015). For instance, the vertical constraint value of 2.0 allows roughly 100 per cent of vertical variation between the constrained parameters. The hydraulic permeability sections are obtained through eq. (5), using directly the inversion parameters σbulk (σ0 for MIC inversions) and $$\sigma _{\rm{max}}^{''}$$ for computation, with σw = 47 mS m − 1. Finally, the depth of investigation (DOI) is computed parameter by parameter following (Fiandaca et al.2015). 3 THE SURVEY 3.1 Survey area The investigated study area (Fig. 2) covers a stretch of Grindsted stream, Denmark. The Grindsted stream is 8–12 m wide and 1–2.5 m deep and flows east to west through the town. It has a discharge of 2000 L/y and drains a sandy aquifer (Balbarini et al.2017; Sonne et al.2017). The groundwater level is very close to the surface near the stream, between 33.5 and 34.0 m a.s.l. (Fig. 2), so the formations are almost completely saturated. A detailed geological description is provided in Section 3.4. At the site, a contaminant plume, consisting of chlorinated ethenes (tetrachlorethene, tricholoroethene, dichloroethene and vinylchloride), BTEX (Benzene, Toluene, Ethylbenzene, Xylenes) and pharmaceutical compounds, is discharging to the stream from the north (Aisopou et al.2015; Rasmussen et al.2016; Balbarini et al.2017; Rønde et al.2017). The investigation was focused around the area north of where the elevated contaminant concentrations were observed in the stream (Sonne et al.2017). Fig. 2 shows the position of the TD SIP profiles (red and blue lines); the multilevel drive point wells used for measuring conductivity in groundwater and hydraulic permeability (green dots); the four wells (B1–B4) used for lithological characterization (black stars), where water conductivity and permeability were also measured. The numbers of the lithological boreholes, which identify the boreholes in the Danish national database (Jupiter), are: 114.2569 (B1); 114.2570 (B2); 114.2507 (B3); 114.1448 (B4). Figure 2. View largeDownload slide Map of the survey area, Grindsted stream, the 2-D TD SIP profiles and wells. Profiles 1–7 are the small-scale 3-D surveys and profiles 8–14 are the large-scale 3-D surveys and profiles 14 and 15 are two additional profiles. The map inserted shows the outline of Denmark, with the location of the Grindsted site marked by the black dot. Figure 2. View largeDownload slide Map of the survey area, Grindsted stream, the 2-D TD SIP profiles and wells. Profiles 1–7 are the small-scale 3-D surveys and profiles 8–14 are the large-scale 3-D surveys and profiles 14 and 15 are two additional profiles. The map inserted shows the outline of Denmark, with the location of the Grindsted site marked by the black dot. 3.2 TD SIP data acquisition and processing TD SIP data were acquired along 16 2-D profiles using a modified Terrameter LS instrument (ABEM). 14 2-D profiles (1–7 and 8–14, see Fig. 2) were collected as part of a large-scale 3-D data survey. The 3-D data acquisition system has been previously demonstrated by Maurya et al. (2017). Each 2-D profile was acquired using a 4 cable layout system. The small-scale profiles (1–7, each 126 m long) were acquired using 64 electrodes with 2 m spacing, where 11 electrodes were connected to the each outer cables (1 and 4) and 21 electrodes were connected to each inner cables (2 and 3). These small-scale profiles were mainly focused around the area where discharge of contaminant is observed in the stream. The long profiles (8–14, each 410 m long) used 63 electrodes each and had 11 electrodes connected to outer cables with 10 m electrode spacing and 21 electrodes connected to the inner cables with 5 m spacing. In these profiles, cable 2 and 3 shares one electrode in the middle of the profile. The interspacing between the big-scale profiles and small-scale profiles was 15 and 7 m, respectively. In addition to the 3-D survey, two additional 2-D profiles (15 and 16 in Fig. 2) of length 600 and 300 m, respectively, were collected before 3-D survey, using the roll-along measurements technique (Dahlin 2001). All the profiles were collected using the gradient array (Dahlin & Zhou 2006) because of the good signal-to-noise ratio of this array (Gazoty et al.2013). Several profiles, both small scale and big scale, cross the meandering of the river (Fig. 2). At these locations, floating electrodes were used. TD SIP data were acquired using 50 per cent duty cycle measurements with current pulse of 4 second on and off time, with full-waveform data recorded at a sampling rate of 3750 Hz. A maximum of 500 mA current was injected with maximum power of 250 watt, however the amount of current was varying (150–500 mA) depending on the contact resistance of the current electrodes. Data were processed and gated into 36 logarithmically spaced gates within the time interval from 1 to 3.9 ms. The Aarhus Workbench software (www.aarhusgeosoftware.dk) was used for manual processing, which mainly involves culling of outliers (individual gates or entire decay) originating from electrode with poor contact. Moreover, electromagnetic induction from the ground or coupling between cables affected the signal at very early time (usually around 1–3 ms) and therefore these data also had to be removed. 3.3 Measurements of water conductivity and hydraulic permeability 3.3.1 Measurements of water conductivity (σw) Water conductivity was measured in groundwater from multi-level drive point wells (10 cm screen lengths and 19 mm inner diameter) and traditional boreholes (1–2 m screen lengths and 50 mm inner diameters). The conductivity was measured at 92 locations, by leading water through a flow-through cell connected to a multimeter (WTW Multi 3420). The conductivity values are measured by the instrument at ambient temperature and then, after temperature correction, recorded at the nominal temperature of 25 °C. Consequently, the conductivity values were then converted to 10 °C according to (Smith 1962), and used for further comparison with the TD SIP data. 3.3.2 Permeability (k) estimation using slug test and grain size analyses To estimate the hydraulic conductivity, falling head slug tests were conducted in the multilevel drive point piezometers and boreholes as the ones described in Section 3.3.1. A pressure transducer, recording the head every 0.1 s, was submerged in the water inside the piezometer/borehole, after which the setup was left to equilibrate. In order to obtain a falling head curve, vacuum was applied, raising the water level ca. 1 m, and then released to let the water level drop (Hinsby et al.1992). The groundwater falling head curves were analysed using Bower and Rice's method (Bouwer & Rice 1976) for screens located in the confined aquifer, Hvorlev's method (Hvorslev 1951) for screens located in the unconfined aquifer, and Springer and Gelhar's method (Springer & Gelhar 1991) for slug tests showing oscillatory responses. Grain size analyses were conducted on soil samples collected every 0.5 m between 0 and 29.5 mbgs at well B3. The grain size distribution curve was determined using sieving, for particle size between 2 and 0.063 mm, and laser diffractometer (Mastersizer Hydro 2000SM), for particle size between 63 and 0.02 μm (Switzer & Pile 2015). The grain size distribution curve was used to estimate the permeability. Many approaches have been suggested for this purpose and the calculated permeability can change orders of magnitudes between the various approaches Devlin (2015). Thus, several methods were applied on each sample (between 2 and 13), depending on the properties of the sample, and the permeability values were computed through the geometric mean. On average, the standard deviation of all methods on all sample is 0.3 decades, indicating that the choice of the qualified methods is not crucial and the proposed approach should give robust estimates of permeability values. 3.4 Geology at the stream site The Grindsted stream is a meandering stream, where sandy sediments are deposited in the channel and banks and peat layers formed on the flood plain. The postglacial stream valley is eroded into a Weichselian outwash plain. The outwash plain is formed by a braided river system located west of the Main stationary Line of the Late Weichselian ice sheet (Houmark-Nielsen 2011). Locally, below the Weichselian outwash plain remains of a Saalian till plain consisting of sand till are found. The sand till overlies a second succession of meltwater deposits, likely of Saalian age (Houmark-Nielsen 2007). The meltwater sand consists mainly of medium grained sand with some beds of silty fine grained sand and gravelly coarse grained sand (Heron et al.1998). The Quaternary deposits of 10–15 m thickness overlie a c. 60 m thick succession of Miocene sediments consisting mainly of mica and quartz sand with some intercalations of clay and lignite beds. In the top part of the Miocene succession, the clay and lignite layers are observed to be up to c. 6 m thick, whereas in the deeper part of the Miocene succession the clay and lignite beds are of c. 1 m thickness. The top and bottom parts of the Miocene succession are dominated by fine grained sand, locally silty and the intermediate parts are coarser grained with medium-coarse grained quartz-rich sand. The sediments belong to the Odderup Formation and are formed in the coastal zone with the clay beds and lignite deposited in lagoonal swamps (Rasmussen et al.2010). The clay and lignite layers can be correlated on a kilometre scale. At c. 80 m below ground surface the top of a clayey succession belonging to the Arnum Formation is situated (Rasmussen et al.2010). A conceptual model of the geological layering in the area (Fig. 3) is generated based on the general understanding of the geological setting (e.g. Heron et al.1998; Rasmussen et al.2010). Approximate depths of the layer boundaries are obtained from borehole lithological logs. Figure 3. View largeDownload slide Conceptual geological model for the Grindsted stream area. Figure 3. View largeDownload slide Conceptual geological model for the Grindsted stream area. 3.5 IP-supported digital 3-D geological model A digital 3-D geological model was constructed for a 0.35 km2 area around Grindsted stream using the modelling software GeoScene3-D (www.i-gis.dk). All available data, including lithological data, water conductivity and hydraulic head observed in boreholes, the imaging results of the IP survey, in terms of σbulk, σ0, $$\sigma _{\rm{max}}^{''}$$ and permeability sections, and a digital terrain model were loaded into the modelling software to be displayed in the 3-D modelling environment. Additionally, historic maps and a geological map were on display in map view. A cognitive modelling approach was followed allowing for incorporation of geological expert knowledge between observational points (Royse 2010). Furthermore, the cognitive modelling process included considerations on methodological limitations of the geophysical information used translating petrophysical parameters into lithological units (Jørgensen et al.2013). With a geological setting consisting of mainly undisturbed layers (Fig. 3), the layer modelling approach was used, where interpretation points initially were placed at locations with data of best quality. When necessary, free interpretation points were added between observational points for constraining surfaces in the interpolation of layer boundaries. The model was divided into a total of 19 layers, 14 Miocene sand and clay/lignite layers, three Quaternary layers consisting of meltwater sand and sand till and two postglacial layers of freshwater sand and peat (Fig. 3). In accordance with the general understanding of the depositional environment, the thin clay and lignite layers are modelled as sub-horizontal continuous layers. The low-permeable clay/lignite layers in the top of Miocene deposits separates an upper aquifer (unconfined) of Quaternary deposits from a lower aquifer (confined) of Miocene deposits. Only three boreholes were available for modelling the top of the lower six Miocene layers (−45 to 0 m a.s.l.) in the deepest part of the model, and additionally five boreholes were present for modelling the top of the next two Miocene layers (0–15 m a.s.l.) at intermediate depths of the model. Within these deeper parts of the 3-D geological model the clay and lignite beds interbedded in the otherwise sandy succession were observed to be c. 1 m thick (Fig. 3), which were beyond the limits of the vertical resolution of the TD SIP method. Thus the layers in this part of the geological model were solely modelled using the few borehole observational points and following the conceptual geological model. In the upper c. 25 m of the geological model (15–40 m a.s.l), all the available data contributed in the modelling process. At the outer parts of the geological model only borehole data were present, whereas in the shallow central and partly contaminated part of the geological model the imaging results of the IP survey (σbulk, σ0 and $$\sigma _{\rm{max}}^{''}$$ and permeability sections) supported the geological modelling in combination with all other available data. However, it was not a trivial task to include the imaging results from IP survey for several reasons: (1) the contamination affected significantly the σ0 profiles, which are classically used for informing geology, although the σbulk, σ''max and permeability sections helped in discriminating contamination and geology; (2) the upper Miocene clay layers of 2–8 m thickness, buried at 12–15 m depth, were at or below the limit of the resolution for the TD SIP method (at least with the acquisition layout used in this study); (3) the petrophysical parameters of the lithological classes, due also to the limited spatial resolution, partly or completely overlapped. This latter point is illustrated in the histogram of Fig. 4, where the permeability values derived from the IP inversions, plotted for three geological units (melt water sand, sand till and clay/lignite), partly overlay each other. The IP-derived permeability values included in the histogram are taken within a horizontal distance of 1.5 times the electrode spacing from the profiles to the boreholes used for the lithological description. Figure 4. View largeDownload slide Histogram of IP-derived permeability values of three geological units: meltwater sand, sand-till and clay/lignite. Figure 4. View largeDownload slide Histogram of IP-derived permeability values of three geological units: meltwater sand, sand-till and clay/lignite. Consequently, if the inverted parameters let to ambiguous lithological interpretations between the borehole data, the geological interpretation then followed the conceptual geological model. 4 RESULTS 4.1 Comparison of BIC and MIC inversions and lithological interpretation Profile 6 was selected for one to one comparison of the MIC and BIC inversion models, because of the presence of a lithological log close to the profile (B3 in Fig. 2), which can help in the result interpretation. This comparison is presented in Fig. 5, where MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). Note that for easy comparison, the total DC conductivity σ0 in the BIC model (Fig. 5b2) was computed using eq. (13). The lithological log is superimposed on the MIC and BIC model parameter sections. DOI is shown for each parameter section by white colour fading, with an upper (more conservative) and lower (less conservative) DOI estimation (Fiandaca et al.2015). Figs 5(f1) and (f2) shows the corresponding data misfit (DC and IP separately) for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo section. It can be noted that both models fits the data equally well, meaning that they resemble equivalent models. However, structural differences can be observed in σ0 and $$\sigma _{\rm{max}}^{''}$$. In the σ0 section computed from the BIC model (Fig. 5b2), the sand-till layer (from depth 7.5 to 10.0 m) and the lignite layer (11.5 to 14 m) are together characterized by a unique, relatively high conductive layer. This also corresponds to a layer with higher imaginary conductivity in $$\sigma _{\rm {max}}^{''}$$ (Fig. 5c2). This is opposed to the MIC model (Fig. 5b1) where the sand-till and lignite layers were not seen as continuous sub-horizontal conductive layers (both in σ0 and $$\sigma _{\rm {max}}^{''}$$). Moreover, the σbulk section of the BIC model shows two relatively homogeneous and distinctive northwest and southwest zones deeper than 10 m. This information is important in discriminating the contribution of conductive lithological material (clay, lignite and sand-till) and water conductivity σw to the DC conductivity. In conclusion we think that inversion results of BIC model parameters best represent the conceptual geological understanding (Fig. 3) and hence our further interpretations are based on this model. The outcome of the joint interpretation of the geophysical results and the lithological and hydrological data is presented in Fig. 5(a1), which represents a slice along profile 6 of the 3-D geological model built using the procedure explained in section 3.5. Figure 5. View largeDownload slide IP inversion results of profile 6, MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). a1 is a slice of the 3-D geological model along profile 6, with lithological log from borehole B3 on top. f1 and f2 show the data misfit for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo-section (blue lines for DC red lines for IP). In f1 and f2, Nite is the number of iterations and χ is the total data misfit. Figure 5. View largeDownload slide IP inversion results of profile 6, MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). a1 is a slice of the 3-D geological model along profile 6, with lithological log from borehole B3 on top. f1 and f2 show the data misfit for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo-section (blue lines for DC red lines for IP). In f1 and f2, Nite is the number of iterations and χ is the total data misfit. 4.2 Mapping of σw, k and lithology in the unconfined aquifer In Fig. 6(a), the correlation between water conductivity (σw) from groundwater and σbulk retrieved from TD SIP inversions is shown. The selection of samples at the stream site for the correlation plot was based on the following three criteria: (1) the value of the model cells closest to the measured sample points were chosen, (2) model cells were included from both 2 and 5 m electrode spacing profiles, however the top 5 m from the 5 m electrode spacing profiles were excluded, considering that the 2 m spacing profiles have better resolution in the near surface and (3) sample points only within the DOI and within an horizontal distance equal to 1.5 times the electrode spacing from profiles were selected. Following the above criteria, a total of 89 sample points out of 92 points were qualified for the comparison. Figure 6. View largeDownload slide (a) Correlation between measured water conductivity,σw and σbulk retrieved from inversion model. Please note that the correlation plot is shown on log–log scale but the fitting of the straight line is performed on a linear scale. (b) Correlation between measured k and IP-derived k values, points filled with light grey colour are the samples taken close to interpreted geological boundaries. The distance from the geological boundaries is the corresponding thickness of the inversion layer. Figure 6. View largeDownload slide (a) Correlation between measured water conductivity,σw and σbulk retrieved from inversion model. Please note that the correlation plot is shown on log–log scale but the fitting of the straight line is performed on a linear scale. (b) Correlation between measured k and IP-derived k values, points filled with light grey colour are the samples taken close to interpreted geological boundaries. The distance from the geological boundaries is the corresponding thickness of the inversion layer. From the correlation plot between σw and σbulk we found the correlation coefficient R2 = 0.31 and the formation factor (inverse of the slope) F = 5.1. The small correlation coefficient is partly due to the limited range of water conductivity measured at the site. In order to get a larger conductivity range, we also plotted in Fig. 6(a) 25 additional points from a recent study by Maurya et al. (2017), which was conducted at a landfill site located 2 km south of the stream in a very similar geological setting. We observed a comparable formation factor (F = 4.7), but a higher correlation coefficient R2 = 0.80 when sample points from landfill are included. Using the value F = 5.1 for the formation factor, water conductivity sections for profile 2, 4, and 6 are shown in Fig. 7. For comparison, the measured water conductivity values are superimposed on the sections. It can be seen that, on average, water conductivities values measured in groundwater agree quite well with the estimates from the TD SIP models. In all the sections, relatively higher water conductivity can be observed in the northwest part of all profiles compared to southwest, which is around the meandering of the river (see Fig. 2). The same criteria as above were adopted for selecting the sample points for investigating the correlation between measured permeability values and those estimated from IP using eq. (5), (Fig. 6b). A total of 94 sample points were qualified for permeability correlation. A fair correlation between measured k values (estimated from both slug test and grain size distribution) and IP-derived k values can be seen in Fig. 6(b). Hydraulic permeability images are produced for profile 2, 4, and 6 using eq. (5), (Fig. 7), with measured k values superimposed on it. Most of the IP-derived k values are within one order of magnitude from the measured values. At few locations (filled with grey colour in Fig. 6b) IP-derived k values are underestimated by more than two orders of magnitude. These samples are located near geological boundaries or in a thin sandy layer interbedded between two low-permeable layers. Considering the vertical smoothness constraint applied in the inversion procedure and the resolution of the TD SIP method, these samples are unlikely to be resolved with the TD SIP method. For example, in profile 6 (Fig. 7d2), the thin sandy layer between sand till and lignite layers (10 to 11.5 m depth interval) is not resolved in the TD SIP inversion. Thus, the TD SIP method underestimates the k in this layer compared to grain size analysis/slug tests. For a quantitative estimation of the prediction quality, the average deviation between the IP-derived k estimates kIP and the measured k values kmeas was computed following the formula $$d = \frac{1}{n}\ \mathop \sum \nolimits_{i\ = \ 1}^n | {{\rm lo{g_{10}}}( {{k_{{\rm{IP}}}}} ) - {\rm lo{g_{10}}}( {{k_{{\rm{meas}}}}} )} |$$ (Weller et al.2015), for all the estimates except the ones near geological boundaries (filled with grey colour in Fig. 6b). The average deviation is d = 0.98 compared to the d = 0.39 value reported by Weller et al. (2015) for laboratory data measured on unconsolidated samples. Except for the discrepancies near geological boundaries, the prediction quality of k values from the IP inversion models is considered satisfactory, both quantitatively and in terms of spatial distribution, and was used to infer the lithology and construct the digital 3-D geological model, which is represented in Fig. 7(a2) along profile 2. The thickening of the sand-till layer from borehole B2 towards the southeast direction comes from the interpretation of the permeability sections, where low-permeability values are seen above 26 m in elevation. Figure 7. View largeDownload slide Zoom-in of the survey map (a1), slice of the geological model along profile 2 (a2), water conductivity (b1–d1) and permeability images (b2–d2) for profiles 2, 4 and 6. Water conductivity sections are derived using σbulkand formation factor of 5.1. Electrical conductivity measured in groundwater and measured permeability values are superimposed on the sections. Only data from wells within 3.0 m from the profile are shown (i.e. 1.5 times the electrode spacing). Note that in profile 6 the continuous distribution of k values is obtained from the grain size analysis. The green dots in panel a1 represent the wells used for water sampling and slug tests/grain size analysis, while the black stars represent the lithological boreholes. The lithological log from borehole B2 is superposed on panel (a2). Figure 7. View largeDownload slide Zoom-in of the survey map (a1), slice of the geological model along profile 2 (a2), water conductivity (b1–d1) and permeability images (b2–d2) for profiles 2, 4 and 6. Water conductivity sections are derived using σbulkand formation factor of 5.1. Electrical conductivity measured in groundwater and measured permeability values are superimposed on the sections. Only data from wells within 3.0 m from the profile are shown (i.e. 1.5 times the electrode spacing). Note that in profile 6 the continuous distribution of k values is obtained from the grain size analysis. The green dots in panel a1 represent the wells used for water sampling and slug tests/grain size analysis, while the black stars represent the lithological boreholes. The lithological log from borehole B2 is superposed on panel (a2). 4.3 Mapping of σw, k and lithology in the confined aquifer The water conductivity and permeability sections of the big-scale profile 14, together with the slice of the 3-D geological model along the profile, are shown in Figs 8(a)–(c), respectively. A relatively higher water conductivity can be seen in the western part of the profile (from 50 to 230 m) below 10 m depth. Similarly, to small-scale profiles, the low-permeability layers from depth of 7 to 14.5 m can be identified as sand till and lignite layers with interbedded meltwater sand. The higher water conductivity is observed in the aquifer below these layers. No contamination was expected in this part of the confined aquifer and confirmation boreholes B1 and B2 were drilled for validating the geophysical results: a good agreement was found with the conductivity of the water samples measured at depths below 20 m, as well as with the permeability estimations at all depths. In particular, the water conductivity measured in the two screens in B1 and B2 at 20 m depth, in the middle of the conductivity anomaly, were 80 mS m − 1 and 60 mS m − 1, well above the 20–30 mS m − 1 background value. For comparison, the values estimated in by the TD SIP inversions in the closest inversion cells, 7 m and 3 m apart, were 120 mS m − 1 and 80 mS m − 1, respectively. Elevated concentrations of inorganic compounds, responsible for the increased water conductivity, were found at B1 and B2 in the deep aquifer. Furthermore, benzene and metabolites of chlorinated solvents were observed in high concentrations (but not as high as in the shallow unconfined aquifer), as well as pharmaceutical contaminants such as sulfonamides and barbiturates. Figure 8. View largeDownload slide Water conductivity, permeability and lithological sections for profile 14. Water sample wells only within 7.5 m from the profile are shown (i.e. within 1.5 times the electrode spacing). Figure 8. View largeDownload slide Water conductivity, permeability and lithological sections for profile 14. Water sample wells only within 7.5 m from the profile are shown (i.e. within 1.5 times the electrode spacing). Finally, Fig. 9 presents pseudo 3-D models of σw and k, obtained by combining all the big-scale sections, and the 3-D geological model. The σw and k models are created by inverse distance interpolations of the 2-D sections to layers and then stacking these layers to construct a 3-D volume. In the σw 3-D model (Fig. 9a), a clear anomaly can be seen localized in the northern part of the survey, whereas in the k 3-D model (Fig. 9b) the low-permeable layer (around 10 − 12.5 to 10 − 14 m2) is more regional and represents the sand-till and lignite/clay layers (Fig. 9c). Figure 9. View largeDownload slide (a) 3-D Geological model, (b) 3-D permeability model and (c) 3-D water electrical conductivity model. Figure 9. View largeDownload slide (a) 3-D Geological model, (b) 3-D permeability model and (c) 3-D water electrical conductivity model. 5 DISCUSSION 5.1 The new approach compared to classical applications of the IP method at contaminated sites For field investigations, risk assessment and adequate remedial design at contaminated sites information on the spatial distribution of the contamination, local geology and hydrogeological properties is required. Aiming expressly at addressing these needs, we developed a new approach for characterizing contaminated sites through TD SIP data. For this purpose, the spatial distributions of two parameters directly usable in the hydrogeological interpretation were derived, that is, the distribution of water conductivity (decoupled from the surface conductivity of the medium) and hydraulic permeability. The water conductivity was used as a proxy for contamination, even though the conductivity response due to contamination is site-/process-specific (e.g. Atekwana & Atekwana 2010) and has not been targeted in this study. On the other hand, the hydraulic permeability has been used as a lithological indicator for mapping the geology (besides its intrinsic value in the hydrogeological characterization). The approach presented in this study differs significantly from the usual applications of IP for characterizing sites impacted by contamination that are generally aimed at the lithological characterization by IP parameters (e.g. Gazoty et al.2012b; Johansson et al.2016) the identification of the source of the contamination (e.g. landfill delineation, Dahlin et al.2010; Gazoty et al.2012a; Wemegah et al.2017), the identification of mobile (pools) or residual NAPLs (e.g. Orozco et al.2012; Johansson et al.2015) or of biogeochemical processes (e.g. Orozco et al.2011; Chen et al.2012). For lithological characterization, parameters related to the magnitude of polarization, such as imaginary conductivity (σ''), phase shift (ϕ) or intrinsic chargeability (m0) are generally used in IP investigations (e.g. Slater & Lesmes 2002b; Gazoty et al.2012b). However, Weller & Slater (2012) have shown that σ'' is dependent on the fluid conductivity, hence lithological interpretation based on σ'' (or ϕ and m0) might be hindered at contaminated sites if significant variability in water conductivity is present. On the contrary, as derived from laboratory results and as shown in this study, permeability estimates of unconsolidated formations depend weakly on water conductivity and hence are a better lithological indicator. Furthermore, while well-established permeability ranges are available for characterizing lithology, it is much more difficult to find appropriate ranges of IP parameters for lithological description in the literature. The proposed approach has an underlying assumption for its applicability: the contamination should have neither IP nor DC signature, except for the effect of water conductivity. This requirement is not always fulfilled, for instance in subsurface settings contaminated with NAPLs (e.g. Orozco et al.2011; Orozco et al.2012; Chen et al.2012; Johansson et al.2015). However, the IP signature is usually significant only where the contaminants are present in concentrations close to the saturation point (e.g. above 1000 mg l−1 for BTEX in Orozco et al.2012), which is generally at or close to the contamination source. Consequently, depending on the contaminant types and concentrations, the distributions of bulk conductivity and hydraulic permeability retrieved by the proposed approach can be inaccurate close to the source area. In contaminant plumes far from the source, the concentrations are usually much lower (typically μg l−1 or few mg l−1 at tens to a few hundreds of meters from the source) and the requirement of a negligible DC/IP signature is entirely fulfilled, as it is also the case at the Grindsted stream site (where for instance the maximum measured BTEX concentration was 1.7 mg l−1). Nevertheless, such plumes may pose a risk to the environment. For instance, in Denmark the limit in groundwater for benzene (a BTEX compound) is 1μg l−1. 5.2 Prediction quality and resolution The use of IP for permeability estimation in the field is not a novelty in itself (e.g. Kemna et al.2004; Hördt et al.2007; Hördt et al.2009; Attwa & Günther 2013). However, to our knowledge the expression suggested by Weller et al. (2015) for permeability estimation (eq. 4), expressly derived for unconsolidated sediments from an extensive set of samples, was used in the field only in the cross-hole survey presented by Binley et al. (2016) and in the logging-while-drilling borehole survey presented by Fiandaca et al. (2017b). Binley et al. (2016) found the field-scale IP method to be suitable for providing estimates of hydraulic permeability in coarse-grained aquifers, but to be somewhat limited for the resolution of small-scale (e.g. lenses versus layers) contrasts in hydraulic permeability variation: low contrast in permeability resulted in relatively low structural resolution based on the distribution of IP parameters. On the contrary, Fiandaca et al. (2017b) found a very good correlation (within on average one order of magnitude) between the IP-derived permeability estimates and those derived using grain size analyses and slug-tests, with similar depth-trends and permeability contrasts. The results presented in this study, obtained only from surface TD SIP measurements, almost replicate the quality of permeability prediction shown in Fiandaca et al. (2017b), except for the lower spatial resolution of the surface imaging that resulted in underestimated predictions close to geological boundaries. Indeed, the spatial resolution of the permeability prediction naturally mimics the spatial resolution of the electrical properties from which the permeability is derived. A way for improving the prediction quality of the TD SIP inversions and to mitigate the resolution issue is to incorporate prior information in the inversion process, both in terms of parameter values and correlations lengths, as done for instance in resistivity inversion by Hermans & Irving (2017). We are currently working on a model parametrization that inverts directly for hydraulic permeability, which allows for a direct integration of the hydrological prior information. The quality of the correlation between the IP-derived hydraulic permeability and the permeability derived from slug tests/grain size analyses is influenced not only by the decrease in resolution with depth of the IP imaging, but more in general by the different support volume of the two estimates. The slug tests, and even more the grain size analyses, provide very local information compared to the surface IP-based permeability estimates. Taking into account also the uncertainty in the petrophysical relationship of eq. (4), this means that a correlation stronger than what we found (Fig. 6b) is most likely not possible. Similar reasoning is valid for the correlation between the water electrical conductivity measured in groundwater and the imaged bulk conductivity (Fig. 6a). The limits in spatial resolution affect also the use of the imaging results for supporting the geological modelling: the geophysical models are sometimes smeared images of the geological models interpreted including all the borehole and geological information. Nevertheless, the imaging of both water electrical conductivity and hydraulic permeability was able to capture the main (hydro) geological units of the site and the areas of higher contamination, and gave a significant input in the site characterization. In particular, this led to the discovery of the contamination in the deep contaminated aquifer. 5.3 Applicability of petrophysical relationships The approach proposed in this study depends on petrophysical relations and it is not applicable when these relations are not valid. In particular, the relationship for permeability estimation (eq. 4) is valid only for unconsolidated, saturated samples. Saturation is also required in the laboratory-derived empirical relation used to isolate the electrolytic bulk conduction from surface conduction (eq. 2). Empirical relations have been suggested in the literature also for the prediction of permeability on consolidated sediments (e.g. Weller et al.2015) and corrections for the dependence on fluid saturation have been proposed for both bulk conductivity (Archie 1942) and surface conductivity (e.g. Vinegar & Waxman 1984; Ulrich & Slater 2004). Nevertheless, the extension of the proposed approach to unsaturated and/or consolidated media is beyond the scope of this study, which deals with saturated unconsolidated sediments. The quality and accuracy of the σbulk estimation through the BIC modelling depends on the value of l used in eq. (2) and on the assumption that the relationship is equally valid for all the investigated sediments. The good correlation found between the imaged bulk conductivity and the water conductivity measured in groundwater in this study (Fig. 6a) is not a proof of the validity of the relationship of eq. (2), but decoupling σbulk from σ'surf led to inversion models more representative of the geological understanding of the site and to reduction of the equivalence problem in the inversions. In fact, the use of the BIC model practically imposes a geometrical constraint between the imaginary conductivity ($$\sigma _{\rm{max}}^{''}$$) section and the total DC conductivity (σ0) section so that a chargeable layer is also conductive. This feature is desirable as long as the petrophysical relation between σ'surf and σ''surf is applicable, because the relation is obeyed by the inversion models by construction, while it might not be fulfilled by inversions carried out for instance with the classic Cole–Cole or the MIC modelling. 5.4 Final remarks The approach presented in this study can contribute to a cost-effective characterization of contaminated sites, reducing the number of drillings needed for the plume delineation and the definition of local geology and hydrogeological properties (also by guiding the drilling campaigns). The learnings on plume, hydraulic permeability field and geology are needed for modelling of contaminant transport, and hence for risk assessment and/or adequate remedial design. For instance, the geological model derived in this study (Fig. 9a), as well as the delineation of the contamination plume (Fig. 9c), will in later studies be used as a framework for building a numerical flow model simulation of contaminant transport towards the stream. 6 CONCLUSION We developed a new approach for characterizing contaminated sites through the imaging of water conductivity σw, hydraulic permeability k, and lithology by means of TD SIP data. We tested the approach at a contaminated stream site using 16 TD SIP profiles and an extensive set of complementary hydrologic and geologic information. For modelling the complex conductivity in the data inversion, a reparametrized Cole–Cole model was developed, namely the bulk and imaginary conductivity (BIC) model, defined in terms of the bulk conductivity σbulk, the MIC $$\sigma _{\rm {max}}^{''}$$, the relaxation time τσ, and the frequency exponent C. In the BIC model the bulk conductivity σbulk is decoupled from the real surface conductivity $$\sigma _{\rm {surf}}^{'}$$ through an empirical, laboratory-derived, petrophysical relation. The subsurface permeability distribution was estimated from the BIC inversion parameters σbulk and σ''max, using a laboratory-derived empirical equation valid for unconsolidated (saturated) sediments without any calibration. The IP-derived permeability values were found to be in good agreement with the estimates obtained using slug tests and grain size analyses: most of the estimates IP-derived k values were on average within one order of magnitude. However, for a few slug tests, made close to geological boundaries, the IP-derived k values were underestimated. This is where the smooth inversion failed to produce sharp boundaries. A positive correlation was observed between water conductivity σw measured in groundwater and the imaged bulk conductivity σbulk, and an average formation factor was estimated to F = 5.1 for converting the imaged values of bulk conductivity into water conductivity. The imaging of permeability and water conductivity allowed for a better discrimination of the clay/lignite lithology from the water conductivity, also due to the decrease in model equivalence of the inversion results, and the geophysical models were actively used for supporting the geological modelling. Furthermore, the direct comparison of the spatial distribution of the imaged water conductivity and the values measured in groundwater, in conjunction with the permeability imaging and the geological interpretation, allowed for the discovery of an unknown increase of σw in the deeper (confined) aquifer, in the northern part of the survey area. Water samples from confirmation wells drilled after the survey showed elevated concentration of inorganic compounds, which are linked to the elevated water conductivity in the confined aquifer. High concentrations of xenobiotic organic contaminants such as benzene, metabolites of chlorinated solvents, sulfonamides, and barbiturates were also observed in the new boreholes. These new findings pave the way for a detailed and inexpensive mapping of bulk/water conductivity, permeability, and lithology at contaminated sites using surface TD SIP measurements, and for cost-effective risk assessment and remedial design. Acknowledgements Support was provided by the research project GEOCON, advancing geological, geophysical and contaminant monitoring technologies for contaminated site investigation (contract 1305-00004B). The funding for GEOCON is provided by Innovation Fund Denmark. REFERENCES Aisopou A., Bjerg P.L., Sonne A.T., Balbarini N., Rosenberg L., Binning P.J., 2015. Dilution and volatilization of groundwater contaminant discharges in streams, J. Contaminant Hydrol. , 172, 71– 83. https://doi.org/10.1016/j.jconhyd.2014.11.004 Google Scholar CrossRef Search ADS   Archie G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. AIME , 146, 54– 62. https://doi.org/10.2118/942054-G Google Scholar CrossRef Search ADS   Atekwana E.A., Atekwana E.A., 2010. Geophysical signatures of microbial activity at hydrocarbon contaminated sites: a review, Surv. Geophys. , 31, 247– 283. https://doi.org/10.1007/s10712-009-9089-8 Google Scholar CrossRef Search ADS   Atekwana E.A., Atekwana E., Legall F.D., Krishnamurthy R.V., 2005. Biodegradation and mineral weathering controls on bulk electrical conductivity in a shallow hydrocarbon contaminated aquifer, J. Contaminant Hydrol. , 80, 149– 167. https://doi.org/10.1016/j.jconhyd.2005.06.009 Google Scholar CrossRef Search ADS   Attwa M., Günther T., 2013. Spectral induced polarization measurements for predicting the hydraulic conductivity in sandy aquifers, Hydrol. Earth Syst. Sci. , 17, 4079– 4094. https://doi.org/10.5194/hess-17-4079-2013 Google Scholar CrossRef Search ADS   Auken E. et al.  , 2015. An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data, Explor. Geophys. , 2015, 223– 235. https://doi.org/10.1071/EG13097 Google Scholar CrossRef Search ADS   Balbarini N., Boon W.M., Nicolajsen E., Nordbotten J.M., Bjerg P.L., Binning P.J., 2017. A 3-D numerical model of the influence of meanders on groundwater discharge to a gaining stream in an unconfined sandy aquifer, J. Hydrol. , 552, 168– 181. https://doi.org/10.1016/j.jhydrol.2017.06.042 Google Scholar CrossRef Search ADS   Barcelona M.J., 1994. Site characterization for subsurface remediation, Ground Water , 32, CONF-9410209. Basu N.B., Rao P., Poyer I.C., Annable M., Hatfield K., 2006. Flux-based assessment at a manufacturing site contaminated with trichloroethylene, J. Contaminant Hydrol. , 86, 105– 127. https://doi.org/10.1016/j.jconhyd.2006.02.011 Google Scholar CrossRef Search ADS   Benson R.C., Yuhr L.B., 2016a. Hydrologic characterization and measurements, in Site Characterization in Karst and Pseudokarst Terraines: Practical Strategies and Technology for Practicing Engineers, Hydrologists and Geologists , pp. 275– 293, eds Benson R.C., Yuhr L.B., Springer Netherlands. Google Scholar CrossRef Search ADS   Benson R.C., Yuhr L.B., 2016b. What is site characterization, in Site Characterization in Karst and Pseudokarst Terraines: Practical Strategies and Technology for Practicing Engineers, Hydrologists and Geologists , pp. 99– 106, eds Benson R.C., Yuhr L.B. Springer Netherlands. Google Scholar CrossRef Search ADS   Binley A., 2015. Tools and Techniques: DC Electrical Methods, in Treatise on Geophysics , 2nd edn, Vol. 11, pp. 233– 259, ed. Schubert G., Elsevier. Google Scholar CrossRef Search ADS   Binley A., Slater L.D., Fukes M., Cassiani G., 2005. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone, Water Resour. Res. , 41, 1– 13. https://doi.org/10.1029/2005WR004202 Google Scholar CrossRef Search ADS   Binley A., Keery J., Slater L., Barrash W., Cardiff M., 2016. The hydrogeologic information in cross-borehole complex conductivity data from an unconsolidated conglomeratic sedimentary aquifer, Geophysics , 81, E409– E421. https://doi.org/10.1190/geo2015-0608.1 Google Scholar CrossRef Search ADS   Bjerg P.L., Tuxen N., Reitzel L.A., Albrechtsen H.-J., Kjeldsen P., 2011. Natural attenuation processes in landfill leachate plumes at three Danish sites, Ground Water , 49, 688– 705. https://doi.org/10.1111/j.1745-6584.2009.00613.x Google Scholar CrossRef Search ADS PubMed  Börner F.D., Schopper J.R., Weller A., 1996. Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements1, Geophys. Prospect. , 44, 583– 601. https://doi.org/10.1111/j.1365-2478.1996.tb00167.x Google Scholar CrossRef Search ADS   Bouwer H., Rice R., 1976. A slug test for determining hydraulic conductivity of unconfined aquifers with completely or partially penetrating wells, Water Resour. Res. , 12, 423– 428. https://doi.org/10.1029/WR012i003p00423 Google Scholar CrossRef Search ADS   Cameron R.E., 1992. Aguide for site and soil description in hazardouswaste site characterization, in Superfund Risk Assessment in Soil Contamination Studies , pp. 3– 17, ed. Hoddinott K.B., ASTM International. Chambers J.E., Kuras O., Meldrum P.I., Ogilvy R.D., Hollands J., 2006. Electrical resistivity tomography applied to geologic, hydrogeologic, and engineering investigations at a former waste-disposal site, Geophysics , 71, B231– B239. https://doi.org/10.1190/1.2360184 Google Scholar CrossRef Search ADS   Chen J., Hubbard S.S., Williams K.H., Orozco A.F., Kemna A., 2012. Estimating the spatiotemporal distribution of geochemical parameters associated with biostimulation using spectral induced polarization data and hierarchical Bayesian models, Water Resour. Res. , 48, 1– 25. https://doi.org/10.1029/2011WR010992 Google Scholar CrossRef Search ADS   Cole K.S., Cole R.H., 1941. Dispersion and absorption in dielectrics I. Alternating current characteristics, J. Chem. Phys. , 9, 341– 351. https://doi.org/10.1063/1.1750906 Google Scholar CrossRef Search ADS   Dahlin T., 2001. The development of DC resistivity imaging techniques, Comput. Geosci. , 27, 1019– 1029. https://doi.org/10.1016/S0098-3004(00)00160-6 Google Scholar CrossRef Search ADS   Dahlin T., Zhou B., 2006. Multiple-gradient array measurements for multichannel 2D resistivity imaging, Near Surf. Geophys. , 4, 113– 123. Dahlin T., Rosqvist H., Leroux V., 2010. Resistivity—IP mapping for landfill applications, First Break , 28, 101– 105. Deryck S.M., Redman J.D., Annan A.P., 1993. Geophysical monitoring of a controlled kerosene spill, in Symposium on the Application of Geophysics to Engineering and Environmental Problems 1993 , pp. 5–19, doi:10.4133/1.2922041. Devlin J.F., 2015. HydrogeoSieveXL: an Excel-based tool to estimate hydraulic conductivity from grain-size analysis, Hydrogeol. J. , 23, 837– 844. https://doi.org/10.1007/s10040-015-1255-0 Google Scholar CrossRef Search ADS   Ellis P.A., Rivett M.O., 2007. Assessing the impact of VOC-contaminated groundwater on surface water at the city scale, J. Contaminant Hydrol. , 91, 107– 127. https://doi.org/10.1016/j.jconhyd.2006.08.015 Google Scholar CrossRef Search ADS   EPA, 1991. Seminar Publication: Site Characterization for Subsurface Remediation , U.S. Environmental Protection Agency. Fiandaca G., Auken E., Gazoty A., Christiansen A.V., 2012. Time-domain-induced polarization: Full-decay forward modeling and 1D laterally constrained inversion of Cole-Cole parameters, Geophysics , 77, E213– E225. https://doi.org/10.1190/geo2011-0217.1 Google Scholar CrossRef Search ADS   Fiandaca G., Ramm J., Binley A., Gazoty A., Christiansen A.V., Auken E., 2013. Resolving spectral information from time domain induced polarization data through 2-D inversion, Geophys. J. Int. , 192, 631– 646. https://doi.org/10.1093/gji/ggs060 Google Scholar CrossRef Search ADS   Fiandaca G., Christiansen A.V., Auken E., 2015. Depth of Investigation for Multi-Parameters Inversions , pp. 666– 670, European Association of Geoscientists & Engineers Publications B.V. (EAGE). Fiandaca G., Madsen L.M., Maurya P.K., 2017a. Re-parameterization of the Cole-Cole Model for Improved Spectral Inversion of Induced Polarization Data, in 23rd European Meeting of Environmental and Engineering Geophysics , doi:10.3997/2214-4609.201702026. Fiandaca G., Maurya P.K., Balbarini N., Hördt A., Christiansen A.V., Foged N., Bjerg P.L., Auken E., 2017b. Hydraulic permeability estimation directly from logging-while-drilling induced polarization data, in AGU-SEG Hydrogeophysics Workshop , Stanford, CA. Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., 2012. Mapping of landfills using time-domain spectral induced polarization data: The Eskelund case study, Near Surf. Geophys. , 10, 575– 586. https://doi.org/10.3997/1873-0604.2012046 Google Scholar CrossRef Search ADS   Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., Pedersen J.K., 2012b. Application of time domain induced polarization to the mapping of lithotypes in a landfill site, Hydrol. Earth Syst. Sci. , 16, 1793– 1804. Google Scholar CrossRef Search ADS   Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., 2013. Data repeatability and acquisition techniques for time-domain spectral induced polarization, Near Surf. Geophys. , 11, 391– 406. https://doi.org/10.3997/1873-0604.2013013 Google Scholar CrossRef Search ADS   Günther T., Martin T., 2016. Spectral two-dimensional inversion of frequency-domain induced polarization data from a mining slag heap, J. Appl. Geophys. , 2016, 10. Hermans T., Irving J., 2017. Facies discrimination with ERT using a probabilistic methodology: effect of sensitivity and regularization, Near Surf. Geophys. , 15( 1), 13– 25. Heron G., Bjerg P.L., Gravesen P., Ludvigsen L., Christensen T.H., 1998. Geology and sediment geochemistry of a landfill leachate contaminated aquifer (Grindsted, Denmark), J. Contaminant Hydrol. , 29, 301– 317. https://doi.org/10.1016/S0169-7722(97)00028-4 Google Scholar CrossRef Search ADS   Hinsby K., Bjerg P.L., Andersen L.J., Skov B., Clausen E.V., 1992. A mini slug test method for determination of a local hydraulic conductivity of an unconfined sandy aquifer, J. Hydrol. , 136, 87– 106. https://doi.org/10.1016/0022-1694(92)90006-H Google Scholar CrossRef Search ADS   Hördt A., Hanstein T., Hönig M., Neubauer F.M., 2006. Efficient spectral IP-modelling in the time domain, J. Appl. Geophys. , 59, 152– 161. https://doi.org/10.1016/j.jappgeo.2005.09.003 Google Scholar CrossRef Search ADS   Hördt A., Blaschek R., Kemna A., Zisser N., 2007. Hydraulic conductivity estimation from induced polarisation data at the field scale — the Krauthausen case history, J. Appl. Geophys. , 62, 33– 46. https://doi.org/10.1016/j.jappgeo.2006.08.001 Google Scholar CrossRef Search ADS   Hördt A., Druiventak A., Blaschek R., Binot F., Kemna A., Kreye P., Zisser N., 2009. Case histories of hydraulic conductivity estimation with induced polarization at the field scale, Near Surf. Geophysics. , 7, 529– 545. https://doi.org/10.3997/1873-0604.2009035 Google Scholar CrossRef Search ADS   Houmark-Nielsen M., 2007. Extent and age of Middle and Late Pleistocene glaciations and periglacial episodes in southern Jylland, Denmark, Bull. Geol. Soc. Denmark , 55, 9– 35. Houmark-Nielsen M., 2011. Pleistocene glaciations in denmark: a closer look at chronology, ice dynamics and landforms, in Developments in Quaternary Sciences , pp. 47– 58, eds Ehlers J., Gibbard P.L., Hughes P.D., Elsevier. Google Scholar CrossRef Search ADS   Hvorslev M., 1951. Time lag and soil permeability in groundwater observations, US Army Corps Eng. Waterways Exp. Sta. Bull. , 36. Hönig M., Tezkan B., 2007. 1D and 2D Cole-Cole-inversion of time-domain induced-polarization data, Geophys. Prospect. , 55, 117– 133. https://doi.org/10.1111/j.1365-2478.2006.00570.x Google Scholar CrossRef Search ADS   Johansson S., Fiandaca G., Dahlin T., 2015. Influence of non-aqueous phase liquid configuration on induced polarization parameters: Conceptual models applied to a time-domain field case study, J. Appl. Geophys. , 123, 295– 309. https://doi.org/10.1016/j.jappgeo.2015.08.010 Google Scholar CrossRef Search ADS   Johansson S., Sparrenbom C., Fiandaca G., Lindskog A., Olsson P.-I., Dahlin T., Rosqvist H., 2017. Investigations of a Cretaceous limestone with spectral induced polarization and scanning electron microscopy, Geophys. J. Int. , 208, 954– 972. https://doi.org/10.1093/gji/ggw432 Google Scholar CrossRef Search ADS   Jørgensen F., Møller R.R., Nebel L., Jensen N., Christiansen A.V., Sandersen P., 2013. A method for cognitive 3D geological voxel modelling of AEM data, Bull. Eng. Geol. Environ. , 72, 421– 432. https://doi.org/10.1007/s10064-013-0487-2 Google Scholar CrossRef Search ADS   Junejo S.A., Zhou Q.Y., Talpur M.A., Debao L., Shaikh S.A., 2015. Imaging of contaminant plumes using ERT in Qinhuai river water and its bed caused by urban effluents at Nanjing, China, Environ. Earth Sci. , 74, 7431– 7440. https://doi.org/10.1007/s12665-015-4728-5 Google Scholar CrossRef Search ADS   Kemna A., Binley A., Slater L., 2004. Crosshole IP imaging for engineering and environmental applications, Geophysics , 69, 97– 107. https://doi.org/10.1190/1.1649379 Google Scholar CrossRef Search ADS   Kemna A. et al.  , 2012. An overview of the spectral induced polarization method for near-surface applications, Near Surf. Geophys. , 10, 453– 468. https://doi.org/10.3997/1873-0604.2012027 Kemna A., Huisman J.A., Zimmermann E., Martin R., Zhao Y., Treichel A., Flores Orozco A., Fechner T., 2014. Broadband electrical impedance tomography for subsurface characterization using improved corrections of electromagnetic coupling and spectral regularization, in Tomography of the Earth's Crust: From Geophysical Sounding to Real-Time Monitoring: GEOTECHNOLOGIEN Science Report No. 21 , pp. 1– 20, eds Weber M., Münch U., Springer International Publishing. Google Scholar CrossRef Search ADS   Lesmes D., Frye K.M., 2001. Influence of pore fluid chemistry on the complex conductivity and induced polarization responses of Berea sandstone, J. geophys. Res. , 160, 4079– 4090. https://doi.org/10.1029/2000JB900392 Google Scholar CrossRef Search ADS   Madsen L.M., Kirkegaard C., Fiandaca G., Christiansen A.V., Auken E., 2016. An analysis of Cole–Cole parameters for IP data using Markov chain Monte Carlo, in IP2016—4th International Workshop on Induced Polarization , Aarhus, Denmark. Madsen L.M., Fiandaca G., Christiansen A.V., Auken E., 2018. Resolution of well-known resistivity equivalences by inclusion of time-domain induced polarization data, Geophysics , 83, E47– E54. Google Scholar CrossRef Search ADS   Maurya P.K., Rønde V.K., Fiandaca G., Balbarini N., Auken E., Bjerg P.L., Christiansen A.V., 2017. Detailed landfill leachate plume mapping using 2D and 3D electrical resistivity tomography—with correlation to ionic strength measured in screens, J. Appl. Geophys. , 138, 1– 8. https://doi.org/10.1016/j.jappgeo.2017.01.019 Google Scholar CrossRef Search ADS   Olsson P.I., Dahlin T., Fiandaca G., Auken E., 2015. Measuring time-domain spectral induced polarization in the on-time:decreasing acquisition time and increasing signal-to-noise ratio, J. Appl. Geophys. , 123, 316– 321. Google Scholar CrossRef Search ADS   Olsson P.-I., Fiandaca G., Larsen J.J., Dahlin T., Auken E., 2016. Doubling the spectrum of time-domain induced polarization by harmonic de-noising, drift correction, spike removal, tapered gating and data uncertainty estimation, Geophys. J. Int. , 207, 774– 784. https://doi.org/10.1093/gji/ggw260 Google Scholar CrossRef Search ADS   Orozco A.F., Williams K.H., Long P.E., Hubbard S.S., Kemna A., 2011. Using complex resistivity imaging to infer biogeochemical processes associated with bioremediation of an uranium-contaminated aquifer, J. geophys. Res. , 116, G03001, doi:10.1029/2010JG001591. Orozco A.F., Kemna A., Oberdörster C., Zschornack L., Leven C., Dietrich P., Weiss H., 2012. Delineation of subsurface hydrocarbon contamination at a former hydrogenation plant using spectral induced polarization imaging, J. Contaminant Hydrol. , 136–137, 131– 144. https://doi.org/10.1016/j.jconhyd.2012.06.001 Google Scholar CrossRef Search ADS   Panagos P., Van Liedekerke M., Yigini Y., Montanarella L., 2013. Contaminated sites in Europe: review of the current situation based on data collected through a European network, J. Environ. Public Health , 2013, 1– 11. https://doi.org/10.1155/2013/158764 Google Scholar CrossRef Search ADS   Pelton W.H., Ward S.H., Hallof P.G., Sill W.R., Nelson P.H., 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP, Geophysics , 43, 588– 609. https://doi.org/10.1190/1.1440839 Google Scholar CrossRef Search ADS   Rasmussen E.S., Dybkjær K., Piasecki S., 2010. Lithostratigraphy of the upper Oligocene - Miocene succession of Denmark, Geol. Surv. Denmark Greenland Bull. , 22, 1– 92. Rasmussen J.J., McKnight U.S., Sonne A.T., Wiberg-Larsen P., Bjerg P.L., 2016. Legacy of a chemical factory site: contaminated groundwater impacts stream macroinvertebrates, Arch. Environ. Contam. Toxicol. , 70, 219– 230. https://doi.org/10.1007/s00244-015-0211-2 Google Scholar CrossRef Search ADS PubMed  Revil A., 2012. Spectral induced polarization of shaly sands: Influence of the electrical double layer, Water Resour. Res. , 48, W02517, doi:10.1029/2011WR011260. https://doi.org/10.1029/2011WR011260 Revil A., Florsch N., 2010. Determination of permeability from spectral induced polarization in granular media, Geophys. J. Int. , 181, 1480– 1498. Rønde V., McKnight U.S., Sonne A.T., Balbarini N., Devlin J.F., Bjerg P.L., 2017. Contaminant mass discharge to streams: Comparing direct groundwater velocity measurements and multi-level groundwater sampling with an in-stream approach, J. Contaminant Hydrol. , 206, 43– 54. https://doi.org/10.1016/j.jconhyd.2017.09.010 Google Scholar CrossRef Search ADS   Royse K.R., 2010. Combining numerical and cognitive 3D modelling approaches in order to determine the structure of the Chalk in the London Basin, Comput. Geosci. , 36, 500– 511. https://doi.org/10.1016/j.cageo.2009.10.001 Google Scholar CrossRef Search ADS   Sims R.C., 1990. Soil remediation techniques at uncontrolled hazardous waste sites, J. Air Waste Manage. Assoc. , 40, 704– 732. https://doi.org/10.1080/10473289.1990.10466716 Google Scholar CrossRef Search ADS PubMed  Slater L., 2007. Near surface electrical characterization of hydraulic conductivity: from petrophysical properties to aquifer geometries-A Review, Surv. Geophys. , 28, 169– 197. https://doi.org/10.1007/s10712-007-9022-y Google Scholar CrossRef Search ADS   Slater L., Lesmes D.P., 2002a. Electrical-hydraulic relationships observed for unconsolidated sediments, Water Resour. Res. , 38, 31-1–31-13. https://doi.org/10.1029/2001WR001075 Slater L.D., Lesmes D., 2002b. IP interpretation in environmental investigations, Geophysics , 67, 77– 88. https://doi.org/10.1190/1.1451353 Google Scholar CrossRef Search ADS   Smith S.H., 1962. Temperature correction in conductivity measurements, Limnol. Oceanogr. , 7, 330– 334. https://doi.org/10.4319/lo.1962.7.3.0330 Google Scholar CrossRef Search ADS   Sonne A.T., McKnight U.S., Rønde V.K., Bjerg P.L., 2017. Assessing the chemical contamination dynamics in a mixed land use stream system, Water Res ., 125, 141– 151. https://doi.org/10.1016/j.watres.2017.08.031 Google Scholar CrossRef Search ADS PubMed  Sparrenbom C.J., Åkesson S., Johansson S., Hagerberg D., Dahlin T., 2017. Investigation of chlorinated solvent pollution with resistivity and induced polarization, Sci. Total Environ. , 575, 767– 778. https://doi.org/10.1016/j.scitotenv.2016.09.117 Google Scholar CrossRef Search ADS PubMed  Springer R., Gelhar L., 1991. Characterization of large-scale aquifer heterogeneity in glacial outwash by analysis of slug tests with oscillatory response, Cape Cod, Massachusetts, US Geol. Surv. Water Res. Invest. Rep. , 91, 36– 40. Switzer A.D., Pile J., 2015. Grain size analysis, in Handbook of Sea-Level Research , eds Shennan I., Long A. J., Horton B. P., John Wiley & Sons, p. 331. Google Scholar CrossRef Search ADS   Sørensen K.I., Larsen F., 1999. Ellog auger drilling: three-in-one method for hydrogeological data collection, Ground Water Monitor Remediation , 19, 97– 101. https://doi.org/10.1111/j.1745-6592.1999.tb00245.x Google Scholar CrossRef Search ADS   Tarasov A., Titov K., 2013. On the use of the Cole–Cole equations in spectral induced polarization, Geophys. J. Int. , 195, 352– 356. https://doi.org/10.1093/gji/ggt251 Google Scholar CrossRef Search ADS   Ulrich C., Slater L.D., 2004. Induced polarization measurements on unsaturated, unconsolidated sands, Geophysics , 69, 762– 771. https://doi.org/10.1190/1.1759462 Google Scholar CrossRef Search ADS   Vinegar H.J., Waxman M.H., 1984. Induced polarization of shaly sands, Geophysics , 49, 1267– 1287. https://doi.org/10.1190/1.1441755 Google Scholar CrossRef Search ADS   Weller A., Slater L., 2012. Salinity dependence of complex conductivity of unconsolidated and consolidated materials: Comparisons with electrical double layer models, Geophysics , 77, D185– D198. https://doi.org/10.1190/geo2012-0030.1 Google Scholar CrossRef Search ADS   Weller A., Breede K., Slater L., Nordsiek S., 2011. Effect of changing water salinity on complex conductivity spectra of sandstones, Geophysics , 76, F315– F327. https://doi.org/10.1190/geo2011-0072.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Nordsiek S., 2013. On the relationship between induced polarization and surface conductivity: implications for petrophysical interpretation of electrical measurements, Geophysics , 78, D315– D325. https://doi.org/10.1190/geo2013-0076.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Binley A., Nordsiek S., Xu S., 2015. Permeability prediction based on induced polarization: insights from measurements on sandstone and unconsolidated samples spanning a wide permeability range, Geophysics , 80, D161– D173. https://doi.org/10.1190/geo2014-0368.1 Google Scholar CrossRef Search ADS   Wemegah D.D., Fiandaca G., Auken E., Menyeh A., Danuor S.K., 2017. Spectral time-domain induced polarisation and magnetic surveying—an efficient tool for characterisation of solid waste deposits in developing countries, Near Surf. Geophys. , 15, 75– 84. Yang C.H., Yu C.Y., Su S.W., 2007. High resistivities associated with a newly formed LNAPL plume imaged by geoelectric techniques—a case study, J. Chin. Inst. Eng. , 30, 53– 62. https://doi.org/10.1080/02533839.2007.9671230 Google Scholar CrossRef Search ADS   Zisser N., Kemna A., Nover G., 2010. Dependence of spectral-induced polarization response of sandstone on temperature and its relevance to permeability estimation, J. geophys. Res. , 115, doi:10.1029/2010JB007526. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Subsurface imaging of water electrical conductivity, hydraulic permeability and lithology at contaminated sites by induced polarization

, Volume 213 (2) – May 1, 2018
16 pages

/lp/ou_press/subsurface-imaging-of-water-electrical-conductivity-hydraulic-qKPpHwHETO
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy018
Publisher site
See Article on Publisher Site

### Abstract

Summary At contaminated sites, knowledge about geology and hydraulic properties of the subsurface and extent of the contamination is needed for assessing the risk and for designing potential site remediation. In this study, we have developed a new approach for characterizing contaminated sites through time-domain spectral induced polarization. The new approach is based on: (1) spectral inversion of the induced polarization data through a reparametrization of the Cole–Cole model, which disentangles the electrolytic bulk conductivity from the surface conductivity for delineating the contamination plume; (2) estimation of hydraulic permeability directly from the inverted parameters using a laboratory-derived empirical equation without any calibration; (3) the use of the geophysical imaging results for supporting the geological modelling and planning of drilling campaigns. The new approach was tested on a data set from the Grindsted stream (Denmark), where contaminated groundwater from a factory site discharges to the stream. Two overlapping areas were covered with seven parallel 2-D profiles each, one large area of 410 m × 90 m (5 m electrode spacing) and one detailed area of 126 m × 42 m (2 m electrode spacing). The geophysical results were complemented and validated by an extensive set of hydrologic and geologic information, including 94 estimates of hydraulic permeability obtained from slug tests and grain size analyses, 89 measurements of water electrical conductivity in groundwater, and four geological logs. On average the IP-derived and measured permeability values agreed within one order of magnitude, except for those close to boundaries between lithological layers (e.g. between sand and clay), where mismatches occurred due to the lack of vertical resolution in the geophysical imaging. An average formation factor was estimated from the correlation between the imaged bulk conductivity values and the water conductivity values measured in groundwater, in order to convert the imaging results from bulk conductivity to water conductivity. The geophysical models were actively used for supporting the geological modelling and the imaging of hydraulic permeability and water conductivity allowed for a better discrimination of the clay/lignite lithology from the pore water conductivity. Furthermore, high water electrical conductivity values were found in a deep confined aquifer, which is separated by a low-permeability clay layer from a shallow aquifer. No contamination was expected in this part of the confined aquifer, and confirmation wells were drilled in the zone of increased water electrical conductivity derived from the geophysical results. Water samples from the new wells showed elevated concentrations of inorganic compounds responsible for the increased water electrical conductivity in the confined aquifer and high concentrations of xenobiotic organic contaminants such as chlorinated ethenes, sulfonamides and barbiturates. Electrical properties, Permeability and porosity, Hydrogeophysics, Inverse theory, Tomography 1 INTRODUCTION Contamination of groundwater and surface water by heavy metals, nutrients or xenobiotic organic compounds in urban areas is a common problem all over the world (Panagos et al.2013). In particular, contaminated sites (e.g. former industrial facilities or dump sites, old landfills, gasoline stations or dry cleaning facilities) can generate groundwater plumes posing a risk to water resources (Basu et al.2006; Ellis & Rivett 2007; Bjerg et al.2011). For site investigations, risk assessment and adequate remedial design, a characterization of the contaminated sites and plumes is needed (Sims 1990; Barcelona 1994). This requires information on the spatial distribution of the contaminant plume, local geology and hydrogeological properties (e.g. permeability, k), which are conventionally estimated from analysis of groundwater samples, of soil samples and of aquifer tests (EPA 1991; Cameron 1992; Benson & Yuhr 2016a). These approaches require many drillings to sufficiently characterize field sites, which is often unfeasible at large sites (Benson & Yuhr 2016b). As an alternative approach, non-invasive geophysical methods can be used in combination with drillings for improved 3-D characterization of contaminant plumes. Among the various geophysical techniques, the direct current (DC) resistivity method has been used for mapping groundwater contamination. The presence of contamination affects the formation resistivity depending on the type and concentration of contaminants. The increase in resistivity is usually associated with the presence of mobile (pools) or residual non-aqueous phase liquids (NAPLs; Deryck et al.1993; Yang et al.2007; Johansson et al.2015). A decrease in the resistivity is typically associated with ionic compounds, which may be linked to contamination and/or to different redox or biodegradation processes (Atekwana et al.2005; Chambers et al.2006; Junejo et al.2015; Maurya et al.2017). The DC resistivity method suffers from well-known limitations, that is, its inability to discriminate between low-resistivity formations (e.g. clay) and high salinity of groundwater. With the IP method, which measures the capacitive nature of the subsurface (Binley 2015), additional information about the surface conductivity is gained (to avoid confusion, here and throughout the manuscript, ‘conductivity’ always refers to electrical conductivity, while ‘permeability’ is used when referring to hydraulic properties). This helps to discriminate lithology-driven resistivity variations from variations driven by pore water. In particular, Weller et al. (2013) identified a strong linear relationship between the real part of surface conductivity and the imaginary conductivity, and discovered how this relationship can be used to improve the estimation of the true formation factor and the groundwater conductivity when an IP measurement is also made. Over the past two decades, the spectral nature of the IP response has increasingly been used as an exploration tool in environmental and hydrogeological investigations (see Kemna et al.2012, for an overview). Spectral IP (SIP) signals can be measured in both time domain (TD) and frequency domain (FD). In TD, SIP measurements are usually combined with DC resistivity measurements using similar field procedures and instrumentation as when measuring DC data alone. Several studies about inversion and processing of TD SIP data have recently been published (Hördt et al.2006; Hönig & Tezkan 2007; Fiandaca et al.2012; Fiandaca et al.2013; Olsson et al.2016; Fiandaca et al.2017a) and successful case studies have been presented where the TD SIP method was applied for mapping lithology (Gazoty et al.2012b; Johansson et al.2016), landfill-waste materials (Gazoty et al.2012a) and contaminated sites (Johansson et al.2015; Sparrenbom et al.2017). In addition to complementing the DC method for lithological discrimination and characterization of contaminated sites, the IP method has been used to estimate hydraulic permeability in the laboratory (e.g. Binley et al.2005; Slater 2007; Revil & Florsch 2010; Zisser et al.2010; Weller et al.2015) and in the field (e.g. Kemna et al.2004; Hördt et al.2007; Hördt et al.2009; Attwa & Günther 2013; Binley et al.2016). In particular, Weller et al. (2015) proposed empirical equations for permeability estimation derived from an extensive set of samples for consolidated and unconsolidated sediments. The applicability of the proposed equation for unconsolidated sediments was verified in the field, in the presence of heterogeneity in both lithology and water chemistry, at the same site presented in this study with the El-log technique (Fiandaca et al.2017b), that is, a borehole technique for acquiring logging-while-drilling DC resistivity, TD spectral IP (SIP) and gamma radiation data (Sørensen & Larsen 1999; Gazoty et al.2012a). The high vertical resolution of the El-log technique matched the lithological variability at the site, minimizing the ambiguity in the result interpretation due to lack of geophysical resolution. A very good correlation within on average one order of magnitude was found between the IP-derived permeability estimates and those derived using grain size analyses and slug-tests, with similar depth-trends and permeability contrasts. These new results on the use of TD SIP data for characterization of hydraulic permeability, the relation between real and imaginary surface conductivity derived in the laboratory, and the recent improvements in processing and inversion of TD SIP data paved the way for the development of a new approach for the characterization of contaminated sites through TD SIP. The new approach is based on: (1) spectral inversion of the TD SIP data through a new reparametrization of the Cole–Cole model (Cole & Cole 1941; Pelton et al.1978), which disentangles the electrolytic bulk conductivity from the surface conductivity for delineating the (inorganic) contamination plume; (2) using the relationship derived from Weller et al. (2015) for estimating the hydraulic permeability in the field directly from the inversion parameters, without any calibration; (3) the use of geophysical imaging results for supporting the geological modelling. The new approach was tested on a TD SIP survey close to the Grindsted stream (Denmark), where the contaminated groundwater from a factory site discharges to the stream (Sonne et al.2017; Rønde et al.2017). The new approach was complemented and validated by comparing the geophysical results with an extensive set of hydrologic and geologic information. In particular, the geophysical results were compared with 94 estimates of hydraulic permeability obtained from slug tests and grain size analyses, 89 measurements of water conductivity in groundwater and four geological logs. 2 METHODOLOGY 2.1 Induced polarization and hydraulic permeability In the induced polarization method the low-frequency capacitive properties of the earth are measured as a frequency-dependent complex resistivity (in FD) or a decay response (in TD) when the medium is excited by a time-varying electric current. In a porous medium free of metallic particles, the electrical conduction takes place through the fluid that fills the interconnected pore spaces (electrolytic bulk conduction) and through the electrical double layer (surface conduction), and these conductivities are usually assumed to add in parallel into the total complex conductivity σ* (e.g. Lesmes & Frye 2001):   \begin{eqnarray} {\sigma ^*} ( f ) &=& {\sigma _{{\rm{bulk}}}}( {{\sigma _{w}}}) + \sigma _{{\rm{surf}}}^* ( {{\sigma _{w},f}}) \nonumber\\ &=& \frac{{{\sigma_{w}}}}{F}\ + \sigma _{{\rm{surf}}}^{'}\left( {{\sigma _{w},f}} \right) + i\sigma _{{\rm{surf}}}^{''}\left( {{\sigma _{w},f}} \right) \end{eqnarray} (1)where the bulk conduction σbulk is expressed by the water conductivity σw divided by the formation factor F and the dependence of the complex surface conductivity $$\sigma _{\rm {surf}}^*$$ over the frequency f and the water conductivity σw is stated expressly. The real and imaginary components of the surface conductivity are not independent: using a database composed of 63 sandstone and unconsolidated sediment samples covering nine independent investigations, Weller et al. (2013) identified a strong linear relationship between $$\sigma _{\rm {surf}}^{'}$$ determined from multisalinity resistivity measurements and $$\sigma _{\rm {surf}}^{''}$$ measured with IP at a frequency of about 1 Hz:   $$\ {\sigma ''_{\rm {surf}}} = \ l \cdot {\sigma '_{\rm {surf}}}$$ (2)with l = 0.042 ± 0.022 (dimension less). Furthermore, Weller et al. (2013) found that the salinity dependency of the real surface conductivity parallels the salinity dependency of the imaginary surface conductivity, which can be expressed as (Weller et al.2011; Weller et al.2015):   $${\sigma ''_{{\rm{surf}}}} ( {{\sigma _w}}) = {\sigma ''_{{\rm{surf}}}}( {{\sigma _f}}) \cdot \frac{1}{{{C_f}}}\sqrt {\frac{{{\sigma _w}}}{{{\sigma _f}}}}$$ (3)where σf represents the salinity of a reference NaCl solution and Cf accounts for the possible differences in ionic species in the reference and actual solution. IP-based permeability prediction methods are based either on the imaginary conductivity $$\sigma _{\rm {surf}}^{''}$$ or on relaxation time τ (a quantity used to represent the characteristic hydraulic length scale e.g. Revil (2012). The premise of estimating k using $$\sigma _{\rm {surf}}^{''}$$ is based on its strong relationship with surface area normalized to the pore volume (Spor), which holds the fundamental basis for derived empirical relationships between k and induced polarization in laboratory studies (Börner et al.1996; Slater & Lesmes 2002a; Weller et al.2015). Recently, Weller et al. (2015) investigated a database consisting of 114 globally collected samples. They avoided the indirect k-estimation through Spor and suggested direct correlations between k and the electrical parameters. For unconsolidated sediments, they proposed the following empirical equation:   $$k\ = \frac{{1.08 \times {{10}^{ - 13}}}}{{{F^{1.12}} \cdot {{\left( {{{\sigma ''}_{{\rm{surf}}}}\left( {{\sigma _f}} \right)} \right)}^{2.27\ }}\ }}.$$ (4)Substituting eq. (3) and using the Archie's law $$({{\sigma _{{\rm{bulk}}}} = \frac{{{\sigma _W}}}{F}})$$, we can re-write eq. (4) as:   $$k\ = \frac{{5.80 \cdot {{10}^{ - 16}}}}{{{C_f}^{2.27}}} \cdot \frac{{{\sigma _{{\rm{bulk}}}}^{1.12}}}{{\sigma {{_{{\rm{surf}}}^{''}}^{2.27\ }}}} \cdot {\sigma _w}^{0.015}$$ (5)where in the last equality the explicit dependence of $$\sigma _{\rm {surf}}^{''}$$ on σw was omitted and the value of σf = 100 mS m − 1 was used. The permeability estimation through eq. (5) depends weakly on water conductivity: a 10-fold and 100-fold variations in water conductivity cause only approximately 3.5 per cent and 7 per cent variations in permeability, respectively. For this reason, in the following the spatial variability of σw is disregarded in permeability computation, and the uniform value σw = 47 mS m − 1 is used, that is, the average value at the site (the range at the site is 11–86 mS m − 1). Weller et al. (2011) suggest Cf = 2 for CaCl2 and Cf = 1 for NaCl; for other ions, no suggestion was made. Considering that a mixture of cations and anions are present in the field-collected water samples with varying molecular concentration, it is difficult to apply an appropriate correction. Therefore, in our k-estimation the value Cf = 1 is used regardless of the chemical composition of the water. 2.2 Parametrization of induced polarization The spectral variation of the complex conductivity is often parametrized, and the Cole–Cole model (Cole & Cole 1941; Pelton et al.1978) is often used for modelling IP data acquired in the field (Fiandaca et al.2013; Kemna et al.2014; Günther & Martin 2016). The Cole–Cole model in its conductivity form is defined as (e.g. Tarasov & Titov 2013)   $${\sigma ^{\rm{*}}} (f) = {\sigma _0}\ \cdot \left[ {1 + \frac{{{m_0}}}{{1 - {m_0}}} \cdot \left( {1 - \frac{1}{{1 + {{\left( {i2\pi f{\tau _\sigma }} \right)}^C}}}} \right)} \right]$$ (6)where σ* is the complex conductivity, σ0 is the DC conductivity, m0 is the intrinsic chargeability, τσ is the relaxation time, C is the frequency exponent, f is the frequency and i is the imaginary unit. Using Monte Carlo Markov Chain analysis, Fiandaca et al. (2017a) have shown that the parameters of the Cole–Cole model are strongly correlated when inverting IP data (in particular m0 and C) and that smaller parameter correlations and better resolution can be achieved (in both FD and TD) when re-parametrizing the Cole–Cole model, replacing the parameter m0 with the maximum imaginary conductivity (MIC) $$\sigma _{{\rm{max}}}^{{\rm{''}}}$$ of the Cole–Cole spectrum (Fig. 1b). Figure 1. View largeDownload slide Spectrum of the BIC model computed using $${\sigma _{{\rm{bulk}}}} = {\rm{\ }}2{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ \ }}\sigma _{{\rm{max}}}^{{\rm{''}}} = {\rm{\ }}0.5{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ }}{\tau _\sigma } = 0.05{\rm{\ s\ and\ }}C = 0.5$$. (a) Real conductivity σ'(black curve), obtained as the sum of the bulk conductivity σbulk (magenta curve) and the surface real conductivity$${\rm{\ }}\sigma _{{\rm{surf}}}^{\rm{'}}$$. The water conductivity value σw with formation factor F = 5 is also shown. (b) Imaginary conductivity σ''. Figure 1. View largeDownload slide Spectrum of the BIC model computed using $${\sigma _{{\rm{bulk}}}} = {\rm{\ }}2{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ \ }}\sigma _{{\rm{max}}}^{{\rm{''}}} = {\rm{\ }}0.5{\rm{\ mS}}\ {{\rm{m}}^{ - 1}},{\rm{\ }}{\tau _\sigma } = 0.05{\rm{\ s\ and\ }}C = 0.5$$. (a) Real conductivity σ'(black curve), obtained as the sum of the bulk conductivity σbulk (magenta curve) and the surface real conductivity$${\rm{\ }}\sigma _{{\rm{surf}}}^{\rm{'}}$$. The water conductivity value σw with formation factor F = 5 is also shown. (b) Imaginary conductivity σ''. The resulting model, named the MIC model, is defined in terms of the parameters mMIC:   $${{\boldsymbol{m}}_{{{\bf MIC}}}} = \{ {{\sigma _0},\sigma _{{\rm{max}}}^{''},\ {\tau _\sigma },\ C} \}$$ (7) In eq. (7) the parameter σ0 represents the total DC conductivity, that is, the sum of the bulk conductivity σbulk and the DC surface conductivity $$\sigma _{{\rm{surf}}}^{'}( {f\ = \ 0} )$$ (see eq. 1), so when imaging the DC conductivity (or, equivalently, DC resistivity), it is difficult to discriminate between low-resistive formations (e.g. clay) and increased salinity of the groundwater. For this reason in this study, we introduce a new Cole–Cole reparametrization, that is, the bulk and (maximum) imaginary conductivity (BIC) model, defined in terms of the parameters mBIC  $${{\boldsymbol{m}}_{{{\bf BIC}}}} = \{ {{\sigma _{{\rm{bulk}}}},\sigma _{{\rm{max}}}^{''},\ {\tau _\sigma },\ C} \}$$ (8)where the assumption is that real and imaginary surface conductivity of eq. (1) are proportional through the proportionality factor l expressed in eq. (2) and that the frequency dependence of the complex conductivity obeys the Cole–Cole model. In Weller et al. (2013), the proportionality factor is found between the DC surface conductivity and the imaginary conductivity at 1 Hz, but the variability of the proportionality with frequency is not discussed. Considering that the imaginary conductivity of the Cole–Cole model reaches a maximum at the frequency f = 1/2πτσ, we decided to enforce in the BIC model the proportionality between the real and imaginary surface conductivity at this frequency. This is a conservative choice: in this way the ratio between the surface imaginary conductivity and real conductivity never exceeds the factor l of eq. (2). On the contrary, enforcing the proportionality at f = 1 Hz would imply a ratio well above l at the peak frequency f = 1/2πτσ for models with τσ ≫ 1. For any given set of BIC model, the corresponding conductivity Cole–Cole parameters can be easy retrieved from the BIC parameters. First, let's define a and b as   $$a\ = \ - {\rm{imag}}\left( {\frac{1}{{1 + {i^C}}}} \right)$$ (9)  $$b\ = \frac{{{m_0}}}{{1 - {m_0}}}\$$ (10) Considering that $${\rm{real}}\ ( {\frac{1}{{1 + {i^C}}}} ) = \ 0.5$$ regardless of the C value, eq. (6) can be written at f = 1/2πτσ as   $${\sigma ^*}( {f\ = 1/2\pi {\tau _\sigma }}) = {\sigma _0}\ \cdot \left[ {1 + b \cdot \left( {0.5 + i \cdot a} \right)} \right]$$ (11) Enforcing the proportionality of eq. (2) at f = 1/2πτσ  and considering that by definition $$\sigma _{{\rm{max}}}^{''} = \ {\rm{imag}}( {{\sigma ^*}( {f\ = 1/2\pi {\tau _\sigma }} )} )$$, we obtain from eq. (1):   $${\sigma _0} \cdot \left[ {1 + b \cdot \left( {0.5 + i \cdot a} \right)} \right]\ = {\sigma _{{\rm{bulk}}}}\ + \sigma _{{\rm{max}}}^{''}/l + \ i \cdot \sigma _{{\rm{max}}}^{''}.$$ (12) The imaginary part of eq. (12) can be solved for σ0:   $${\sigma _0} = \frac{\sigma _{\rm{max}}^{''}} {a \cdot b}.$$ (13) Substituting eq. (13) into eq. (12), it is possible to solve for b from the real part of the resulting equation:   $$b = \frac{l \cdot \sigma _{\rm{max}}^{''}/a} {{\sigma _{\rm{max}}^{''}} + l \cdot {\sigma _{\rm{bulk}}} - 0.5 \cdot l \cdot \sigma _{\rm{max}}^{''}/a}$$ (14) Finally, m0 can be retrieved as a function of b from eq. (10):   $${m_0} = \frac{b}{{1 + b}}\$$ (15) For example, Fig. 1 shows the spectrum of the BIC model defined by the parameter values $$\{ \sigma _{\rm{bulk}} = 2\ {\rm{mS}}\ {\rm{m}}^{ - 1},\ \sigma _{\rm{max}}^{''} = 0.5\ {\rm{mS}}\ {\rm{m}}^{ - 1},\ \tau _\sigma = 0.05\ {\rm{s}}, C = 0.5 \}$$. The σ0 and m0 values of the corresponding Cole–Cole model are σ0 = 12.7 mS m − 1 and m0 = 160 mV V − 1, and the DC conductivity is almost completely dominated by the surface conductivity. 2.3 Time-domain spectral induced polarization TD spectral IP is an extension of the DC resistivity method in which the capacity of the ground is measured using square current pulses. A TD SIP signal can either be measured as rising of the potential during the current on time or decaying of potential after the current is turned off. The measurement using the former approach is referred as 100 per cent duty cycle measurement and the latter as 50 per cent duty cycle measurement (Olsson et al.2015). The possibility of extracting the spectral information contained in the TD signal depends considerably on the number of decades acquired (Madsen et al.2016, 2018), and nowadays commercial instruments exist that sample the full-waveform potential and current signals at high sampling rate (typically in the kHz range) for the entire measurement time, allowing for advanced signal processing. The recorded full-waveform data often contain harmonic noise from the power distribution grid, spikes from, for example, animal fences and self-potential background drift. Harmonic noise makes it impossible to measure earlier than 20 ms (in a 50 Hz environment) and self-potential drift distorts the signal at later times: this limits the extraction of spectral information from the signal. In this study full-waveform data sampled at 3750 Hz were acquired and processed following Olsson et al. (2016) for removal of harmonic noise, spikes, self-potential drift and tapered windowing, in order to retrieve unbiased TD SIP data from a few milliseconds after current switch. Apparent DC resistivity values and full-decay gated curves (i.e. apparent chargeability values with time) represent the data space for the inversion of TD SIP data. 2.4 2-D inversion and imaging of hydraulic permeability The inversion procedure is carried out using the code by Auken et al. (2015), where the apparent resistivity values and the IP full-decays are inverted simultaneously in 2-D following Fiandaca et al. (2013), with an accurate description of transmitter waveform and the receiver transfer function for an unbiased estimation of the spectral parameters (Fiandaca et al.2012). Along the 2-D section, the model space is defined cell by cell in terms of a parametrization of the FD complex conductivity (The MIC and BIC parametrizations are the ones used in this study). The same vertical model discretization (with an increased number of layers in the big-scale profiles), model constraints and starting model were used for all the inversions, as well as the same stopping criterion for the iterative inversion (i.e. a relative variation of the objective function between consecutive iterations below 2 per cent). Furthermore, the same model for the data standard deviations (STD) has been adopted, that is, 1 per cent on apparent resistivity and 10 per cent on chargeability, plus a noise floor of 0.1 mV (Olsson et al.2015). L2 smoothness vertical/horizontal constraints were applied, with values 2.0 and 1.2, respectively. The constraint values represent the relative vertical/lateral variation of the parameters that weights the roughness misfit in the objective function (Auken et al.2015). For instance, the vertical constraint value of 2.0 allows roughly 100 per cent of vertical variation between the constrained parameters. The hydraulic permeability sections are obtained through eq. (5), using directly the inversion parameters σbulk (σ0 for MIC inversions) and $$\sigma _{\rm{max}}^{''}$$ for computation, with σw = 47 mS m − 1. Finally, the depth of investigation (DOI) is computed parameter by parameter following (Fiandaca et al.2015). 3 THE SURVEY 3.1 Survey area The investigated study area (Fig. 2) covers a stretch of Grindsted stream, Denmark. The Grindsted stream is 8–12 m wide and 1–2.5 m deep and flows east to west through the town. It has a discharge of 2000 L/y and drains a sandy aquifer (Balbarini et al.2017; Sonne et al.2017). The groundwater level is very close to the surface near the stream, between 33.5 and 34.0 m a.s.l. (Fig. 2), so the formations are almost completely saturated. A detailed geological description is provided in Section 3.4. At the site, a contaminant plume, consisting of chlorinated ethenes (tetrachlorethene, tricholoroethene, dichloroethene and vinylchloride), BTEX (Benzene, Toluene, Ethylbenzene, Xylenes) and pharmaceutical compounds, is discharging to the stream from the north (Aisopou et al.2015; Rasmussen et al.2016; Balbarini et al.2017; Rønde et al.2017). The investigation was focused around the area north of where the elevated contaminant concentrations were observed in the stream (Sonne et al.2017). Fig. 2 shows the position of the TD SIP profiles (red and blue lines); the multilevel drive point wells used for measuring conductivity in groundwater and hydraulic permeability (green dots); the four wells (B1–B4) used for lithological characterization (black stars), where water conductivity and permeability were also measured. The numbers of the lithological boreholes, which identify the boreholes in the Danish national database (Jupiter), are: 114.2569 (B1); 114.2570 (B2); 114.2507 (B3); 114.1448 (B4). Figure 2. View largeDownload slide Map of the survey area, Grindsted stream, the 2-D TD SIP profiles and wells. Profiles 1–7 are the small-scale 3-D surveys and profiles 8–14 are the large-scale 3-D surveys and profiles 14 and 15 are two additional profiles. The map inserted shows the outline of Denmark, with the location of the Grindsted site marked by the black dot. Figure 2. View largeDownload slide Map of the survey area, Grindsted stream, the 2-D TD SIP profiles and wells. Profiles 1–7 are the small-scale 3-D surveys and profiles 8–14 are the large-scale 3-D surveys and profiles 14 and 15 are two additional profiles. The map inserted shows the outline of Denmark, with the location of the Grindsted site marked by the black dot. 3.2 TD SIP data acquisition and processing TD SIP data were acquired along 16 2-D profiles using a modified Terrameter LS instrument (ABEM). 14 2-D profiles (1–7 and 8–14, see Fig. 2) were collected as part of a large-scale 3-D data survey. The 3-D data acquisition system has been previously demonstrated by Maurya et al. (2017). Each 2-D profile was acquired using a 4 cable layout system. The small-scale profiles (1–7, each 126 m long) were acquired using 64 electrodes with 2 m spacing, where 11 electrodes were connected to the each outer cables (1 and 4) and 21 electrodes were connected to each inner cables (2 and 3). These small-scale profiles were mainly focused around the area where discharge of contaminant is observed in the stream. The long profiles (8–14, each 410 m long) used 63 electrodes each and had 11 electrodes connected to outer cables with 10 m electrode spacing and 21 electrodes connected to the inner cables with 5 m spacing. In these profiles, cable 2 and 3 shares one electrode in the middle of the profile. The interspacing between the big-scale profiles and small-scale profiles was 15 and 7 m, respectively. In addition to the 3-D survey, two additional 2-D profiles (15 and 16 in Fig. 2) of length 600 and 300 m, respectively, were collected before 3-D survey, using the roll-along measurements technique (Dahlin 2001). All the profiles were collected using the gradient array (Dahlin & Zhou 2006) because of the good signal-to-noise ratio of this array (Gazoty et al.2013). Several profiles, both small scale and big scale, cross the meandering of the river (Fig. 2). At these locations, floating electrodes were used. TD SIP data were acquired using 50 per cent duty cycle measurements with current pulse of 4 second on and off time, with full-waveform data recorded at a sampling rate of 3750 Hz. A maximum of 500 mA current was injected with maximum power of 250 watt, however the amount of current was varying (150–500 mA) depending on the contact resistance of the current electrodes. Data were processed and gated into 36 logarithmically spaced gates within the time interval from 1 to 3.9 ms. The Aarhus Workbench software (www.aarhusgeosoftware.dk) was used for manual processing, which mainly involves culling of outliers (individual gates or entire decay) originating from electrode with poor contact. Moreover, electromagnetic induction from the ground or coupling between cables affected the signal at very early time (usually around 1–3 ms) and therefore these data also had to be removed. 3.3 Measurements of water conductivity and hydraulic permeability 3.3.1 Measurements of water conductivity (σw) Water conductivity was measured in groundwater from multi-level drive point wells (10 cm screen lengths and 19 mm inner diameter) and traditional boreholes (1–2 m screen lengths and 50 mm inner diameters). The conductivity was measured at 92 locations, by leading water through a flow-through cell connected to a multimeter (WTW Multi 3420). The conductivity values are measured by the instrument at ambient temperature and then, after temperature correction, recorded at the nominal temperature of 25 °C. Consequently, the conductivity values were then converted to 10 °C according to (Smith 1962), and used for further comparison with the TD SIP data. 3.3.2 Permeability (k) estimation using slug test and grain size analyses To estimate the hydraulic conductivity, falling head slug tests were conducted in the multilevel drive point piezometers and boreholes as the ones described in Section 3.3.1. A pressure transducer, recording the head every 0.1 s, was submerged in the water inside the piezometer/borehole, after which the setup was left to equilibrate. In order to obtain a falling head curve, vacuum was applied, raising the water level ca. 1 m, and then released to let the water level drop (Hinsby et al.1992). The groundwater falling head curves were analysed using Bower and Rice's method (Bouwer & Rice 1976) for screens located in the confined aquifer, Hvorlev's method (Hvorslev 1951) for screens located in the unconfined aquifer, and Springer and Gelhar's method (Springer & Gelhar 1991) for slug tests showing oscillatory responses. Grain size analyses were conducted on soil samples collected every 0.5 m between 0 and 29.5 mbgs at well B3. The grain size distribution curve was determined using sieving, for particle size between 2 and 0.063 mm, and laser diffractometer (Mastersizer Hydro 2000SM), for particle size between 63 and 0.02 μm (Switzer & Pile 2015). The grain size distribution curve was used to estimate the permeability. Many approaches have been suggested for this purpose and the calculated permeability can change orders of magnitudes between the various approaches Devlin (2015). Thus, several methods were applied on each sample (between 2 and 13), depending on the properties of the sample, and the permeability values were computed through the geometric mean. On average, the standard deviation of all methods on all sample is 0.3 decades, indicating that the choice of the qualified methods is not crucial and the proposed approach should give robust estimates of permeability values. 3.4 Geology at the stream site The Grindsted stream is a meandering stream, where sandy sediments are deposited in the channel and banks and peat layers formed on the flood plain. The postglacial stream valley is eroded into a Weichselian outwash plain. The outwash plain is formed by a braided river system located west of the Main stationary Line of the Late Weichselian ice sheet (Houmark-Nielsen 2011). Locally, below the Weichselian outwash plain remains of a Saalian till plain consisting of sand till are found. The sand till overlies a second succession of meltwater deposits, likely of Saalian age (Houmark-Nielsen 2007). The meltwater sand consists mainly of medium grained sand with some beds of silty fine grained sand and gravelly coarse grained sand (Heron et al.1998). The Quaternary deposits of 10–15 m thickness overlie a c. 60 m thick succession of Miocene sediments consisting mainly of mica and quartz sand with some intercalations of clay and lignite beds. In the top part of the Miocene succession, the clay and lignite layers are observed to be up to c. 6 m thick, whereas in the deeper part of the Miocene succession the clay and lignite beds are of c. 1 m thickness. The top and bottom parts of the Miocene succession are dominated by fine grained sand, locally silty and the intermediate parts are coarser grained with medium-coarse grained quartz-rich sand. The sediments belong to the Odderup Formation and are formed in the coastal zone with the clay beds and lignite deposited in lagoonal swamps (Rasmussen et al.2010). The clay and lignite layers can be correlated on a kilometre scale. At c. 80 m below ground surface the top of a clayey succession belonging to the Arnum Formation is situated (Rasmussen et al.2010). A conceptual model of the geological layering in the area (Fig. 3) is generated based on the general understanding of the geological setting (e.g. Heron et al.1998; Rasmussen et al.2010). Approximate depths of the layer boundaries are obtained from borehole lithological logs. Figure 3. View largeDownload slide Conceptual geological model for the Grindsted stream area. Figure 3. View largeDownload slide Conceptual geological model for the Grindsted stream area. 3.5 IP-supported digital 3-D geological model A digital 3-D geological model was constructed for a 0.35 km2 area around Grindsted stream using the modelling software GeoScene3-D (www.i-gis.dk). All available data, including lithological data, water conductivity and hydraulic head observed in boreholes, the imaging results of the IP survey, in terms of σbulk, σ0, $$\sigma _{\rm{max}}^{''}$$ and permeability sections, and a digital terrain model were loaded into the modelling software to be displayed in the 3-D modelling environment. Additionally, historic maps and a geological map were on display in map view. A cognitive modelling approach was followed allowing for incorporation of geological expert knowledge between observational points (Royse 2010). Furthermore, the cognitive modelling process included considerations on methodological limitations of the geophysical information used translating petrophysical parameters into lithological units (Jørgensen et al.2013). With a geological setting consisting of mainly undisturbed layers (Fig. 3), the layer modelling approach was used, where interpretation points initially were placed at locations with data of best quality. When necessary, free interpretation points were added between observational points for constraining surfaces in the interpolation of layer boundaries. The model was divided into a total of 19 layers, 14 Miocene sand and clay/lignite layers, three Quaternary layers consisting of meltwater sand and sand till and two postglacial layers of freshwater sand and peat (Fig. 3). In accordance with the general understanding of the depositional environment, the thin clay and lignite layers are modelled as sub-horizontal continuous layers. The low-permeable clay/lignite layers in the top of Miocene deposits separates an upper aquifer (unconfined) of Quaternary deposits from a lower aquifer (confined) of Miocene deposits. Only three boreholes were available for modelling the top of the lower six Miocene layers (−45 to 0 m a.s.l.) in the deepest part of the model, and additionally five boreholes were present for modelling the top of the next two Miocene layers (0–15 m a.s.l.) at intermediate depths of the model. Within these deeper parts of the 3-D geological model the clay and lignite beds interbedded in the otherwise sandy succession were observed to be c. 1 m thick (Fig. 3), which were beyond the limits of the vertical resolution of the TD SIP method. Thus the layers in this part of the geological model were solely modelled using the few borehole observational points and following the conceptual geological model. In the upper c. 25 m of the geological model (15–40 m a.s.l), all the available data contributed in the modelling process. At the outer parts of the geological model only borehole data were present, whereas in the shallow central and partly contaminated part of the geological model the imaging results of the IP survey (σbulk, σ0 and $$\sigma _{\rm{max}}^{''}$$ and permeability sections) supported the geological modelling in combination with all other available data. However, it was not a trivial task to include the imaging results from IP survey for several reasons: (1) the contamination affected significantly the σ0 profiles, which are classically used for informing geology, although the σbulk, σ''max and permeability sections helped in discriminating contamination and geology; (2) the upper Miocene clay layers of 2–8 m thickness, buried at 12–15 m depth, were at or below the limit of the resolution for the TD SIP method (at least with the acquisition layout used in this study); (3) the petrophysical parameters of the lithological classes, due also to the limited spatial resolution, partly or completely overlapped. This latter point is illustrated in the histogram of Fig. 4, where the permeability values derived from the IP inversions, plotted for three geological units (melt water sand, sand till and clay/lignite), partly overlay each other. The IP-derived permeability values included in the histogram are taken within a horizontal distance of 1.5 times the electrode spacing from the profiles to the boreholes used for the lithological description. Figure 4. View largeDownload slide Histogram of IP-derived permeability values of three geological units: meltwater sand, sand-till and clay/lignite. Figure 4. View largeDownload slide Histogram of IP-derived permeability values of three geological units: meltwater sand, sand-till and clay/lignite. Consequently, if the inverted parameters let to ambiguous lithological interpretations between the borehole data, the geological interpretation then followed the conceptual geological model. 4 RESULTS 4.1 Comparison of BIC and MIC inversions and lithological interpretation Profile 6 was selected for one to one comparison of the MIC and BIC inversion models, because of the presence of a lithological log close to the profile (B3 in Fig. 2), which can help in the result interpretation. This comparison is presented in Fig. 5, where MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). Note that for easy comparison, the total DC conductivity σ0 in the BIC model (Fig. 5b2) was computed using eq. (13). The lithological log is superimposed on the MIC and BIC model parameter sections. DOI is shown for each parameter section by white colour fading, with an upper (more conservative) and lower (less conservative) DOI estimation (Fiandaca et al.2015). Figs 5(f1) and (f2) shows the corresponding data misfit (DC and IP separately) for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo section. It can be noted that both models fits the data equally well, meaning that they resemble equivalent models. However, structural differences can be observed in σ0 and $$\sigma _{\rm{max}}^{''}$$. In the σ0 section computed from the BIC model (Fig. 5b2), the sand-till layer (from depth 7.5 to 10.0 m) and the lignite layer (11.5 to 14 m) are together characterized by a unique, relatively high conductive layer. This also corresponds to a layer with higher imaginary conductivity in $$\sigma _{\rm {max}}^{''}$$ (Fig. 5c2). This is opposed to the MIC model (Fig. 5b1) where the sand-till and lignite layers were not seen as continuous sub-horizontal conductive layers (both in σ0 and $$\sigma _{\rm {max}}^{''}$$). Moreover, the σbulk section of the BIC model shows two relatively homogeneous and distinctive northwest and southwest zones deeper than 10 m. This information is important in discriminating the contribution of conductive lithological material (clay, lignite and sand-till) and water conductivity σw to the DC conductivity. In conclusion we think that inversion results of BIC model parameters best represent the conceptual geological understanding (Fig. 3) and hence our further interpretations are based on this model. The outcome of the joint interpretation of the geophysical results and the lithological and hydrological data is presented in Fig. 5(a1), which represents a slice along profile 6 of the 3-D geological model built using the procedure explained in section 3.5. Figure 5. View largeDownload slide IP inversion results of profile 6, MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). a1 is a slice of the 3-D geological model along profile 6, with lithological log from borehole B3 on top. f1 and f2 show the data misfit for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo-section (blue lines for DC red lines for IP). In f1 and f2, Nite is the number of iterations and χ is the total data misfit. Figure 5. View largeDownload slide IP inversion results of profile 6, MIC model parameters are shown in left panel (b1–e1) and BIC model parameters are shown in right panel (a2–e2). a1 is a slice of the 3-D geological model along profile 6, with lithological log from borehole B3 on top. f1 and f2 show the data misfit for the MIC and BIC inversions, averaged vertically (and over all gates for the IP misfit) along the pseudo-section (blue lines for DC red lines for IP). In f1 and f2, Nite is the number of iterations and χ is the total data misfit. 4.2 Mapping of σw, k and lithology in the unconfined aquifer In Fig. 6(a), the correlation between water conductivity (σw) from groundwater and σbulk retrieved from TD SIP inversions is shown. The selection of samples at the stream site for the correlation plot was based on the following three criteria: (1) the value of the model cells closest to the measured sample points were chosen, (2) model cells were included from both 2 and 5 m electrode spacing profiles, however the top 5 m from the 5 m electrode spacing profiles were excluded, considering that the 2 m spacing profiles have better resolution in the near surface and (3) sample points only within the DOI and within an horizontal distance equal to 1.5 times the electrode spacing from profiles were selected. Following the above criteria, a total of 89 sample points out of 92 points were qualified for the comparison. Figure 6. View largeDownload slide (a) Correlation between measured water conductivity,σw and σbulk retrieved from inversion model. Please note that the correlation plot is shown on log–log scale but the fitting of the straight line is performed on a linear scale. (b) Correlation between measured k and IP-derived k values, points filled with light grey colour are the samples taken close to interpreted geological boundaries. The distance from the geological boundaries is the corresponding thickness of the inversion layer. Figure 6. View largeDownload slide (a) Correlation between measured water conductivity,σw and σbulk retrieved from inversion model. Please note that the correlation plot is shown on log–log scale but the fitting of the straight line is performed on a linear scale. (b) Correlation between measured k and IP-derived k values, points filled with light grey colour are the samples taken close to interpreted geological boundaries. The distance from the geological boundaries is the corresponding thickness of the inversion layer. From the correlation plot between σw and σbulk we found the correlation coefficient R2 = 0.31 and the formation factor (inverse of the slope) F = 5.1. The small correlation coefficient is partly due to the limited range of water conductivity measured at the site. In order to get a larger conductivity range, we also plotted in Fig. 6(a) 25 additional points from a recent study by Maurya et al. (2017), which was conducted at a landfill site located 2 km south of the stream in a very similar geological setting. We observed a comparable formation factor (F = 4.7), but a higher correlation coefficient R2 = 0.80 when sample points from landfill are included. Using the value F = 5.1 for the formation factor, water conductivity sections for profile 2, 4, and 6 are shown in Fig. 7. For comparison, the measured water conductivity values are superimposed on the sections. It can be seen that, on average, water conductivities values measured in groundwater agree quite well with the estimates from the TD SIP models. In all the sections, relatively higher water conductivity can be observed in the northwest part of all profiles compared to southwest, which is around the meandering of the river (see Fig. 2). The same criteria as above were adopted for selecting the sample points for investigating the correlation between measured permeability values and those estimated from IP using eq. (5), (Fig. 6b). A total of 94 sample points were qualified for permeability correlation. A fair correlation between measured k values (estimated from both slug test and grain size distribution) and IP-derived k values can be seen in Fig. 6(b). Hydraulic permeability images are produced for profile 2, 4, and 6 using eq. (5), (Fig. 7), with measured k values superimposed on it. Most of the IP-derived k values are within one order of magnitude from the measured values. At few locations (filled with grey colour in Fig. 6b) IP-derived k values are underestimated by more than two orders of magnitude. These samples are located near geological boundaries or in a thin sandy layer interbedded between two low-permeable layers. Considering the vertical smoothness constraint applied in the inversion procedure and the resolution of the TD SIP method, these samples are unlikely to be resolved with the TD SIP method. For example, in profile 6 (Fig. 7d2), the thin sandy layer between sand till and lignite layers (10 to 11.5 m depth interval) is not resolved in the TD SIP inversion. Thus, the TD SIP method underestimates the k in this layer compared to grain size analysis/slug tests. For a quantitative estimation of the prediction quality, the average deviation between the IP-derived k estimates kIP and the measured k values kmeas was computed following the formula $$d = \frac{1}{n}\ \mathop \sum \nolimits_{i\ = \ 1}^n | {{\rm lo{g_{10}}}( {{k_{{\rm{IP}}}}} ) - {\rm lo{g_{10}}}( {{k_{{\rm{meas}}}}} )} |$$ (Weller et al.2015), for all the estimates except the ones near geological boundaries (filled with grey colour in Fig. 6b). The average deviation is d = 0.98 compared to the d = 0.39 value reported by Weller et al. (2015) for laboratory data measured on unconsolidated samples. Except for the discrepancies near geological boundaries, the prediction quality of k values from the IP inversion models is considered satisfactory, both quantitatively and in terms of spatial distribution, and was used to infer the lithology and construct the digital 3-D geological model, which is represented in Fig. 7(a2) along profile 2. The thickening of the sand-till layer from borehole B2 towards the southeast direction comes from the interpretation of the permeability sections, where low-permeability values are seen above 26 m in elevation. Figure 7. View largeDownload slide Zoom-in of the survey map (a1), slice of the geological model along profile 2 (a2), water conductivity (b1–d1) and permeability images (b2–d2) for profiles 2, 4 and 6. Water conductivity sections are derived using σbulkand formation factor of 5.1. Electrical conductivity measured in groundwater and measured permeability values are superimposed on the sections. Only data from wells within 3.0 m from the profile are shown (i.e. 1.5 times the electrode spacing). Note that in profile 6 the continuous distribution of k values is obtained from the grain size analysis. The green dots in panel a1 represent the wells used for water sampling and slug tests/grain size analysis, while the black stars represent the lithological boreholes. The lithological log from borehole B2 is superposed on panel (a2). Figure 7. View largeDownload slide Zoom-in of the survey map (a1), slice of the geological model along profile 2 (a2), water conductivity (b1–d1) and permeability images (b2–d2) for profiles 2, 4 and 6. Water conductivity sections are derived using σbulkand formation factor of 5.1. Electrical conductivity measured in groundwater and measured permeability values are superimposed on the sections. Only data from wells within 3.0 m from the profile are shown (i.e. 1.5 times the electrode spacing). Note that in profile 6 the continuous distribution of k values is obtained from the grain size analysis. The green dots in panel a1 represent the wells used for water sampling and slug tests/grain size analysis, while the black stars represent the lithological boreholes. The lithological log from borehole B2 is superposed on panel (a2). 4.3 Mapping of σw, k and lithology in the confined aquifer The water conductivity and permeability sections of the big-scale profile 14, together with the slice of the 3-D geological model along the profile, are shown in Figs 8(a)–(c), respectively. A relatively higher water conductivity can be seen in the western part of the profile (from 50 to 230 m) below 10 m depth. Similarly, to small-scale profiles, the low-permeability layers from depth of 7 to 14.5 m can be identified as sand till and lignite layers with interbedded meltwater sand. The higher water conductivity is observed in the aquifer below these layers. No contamination was expected in this part of the confined aquifer and confirmation boreholes B1 and B2 were drilled for validating the geophysical results: a good agreement was found with the conductivity of the water samples measured at depths below 20 m, as well as with the permeability estimations at all depths. In particular, the water conductivity measured in the two screens in B1 and B2 at 20 m depth, in the middle of the conductivity anomaly, were 80 mS m − 1 and 60 mS m − 1, well above the 20–30 mS m − 1 background value. For comparison, the values estimated in by the TD SIP inversions in the closest inversion cells, 7 m and 3 m apart, were 120 mS m − 1 and 80 mS m − 1, respectively. Elevated concentrations of inorganic compounds, responsible for the increased water conductivity, were found at B1 and B2 in the deep aquifer. Furthermore, benzene and metabolites of chlorinated solvents were observed in high concentrations (but not as high as in the shallow unconfined aquifer), as well as pharmaceutical contaminants such as sulfonamides and barbiturates. Figure 8. View largeDownload slide Water conductivity, permeability and lithological sections for profile 14. Water sample wells only within 7.5 m from the profile are shown (i.e. within 1.5 times the electrode spacing). Figure 8. View largeDownload slide Water conductivity, permeability and lithological sections for profile 14. Water sample wells only within 7.5 m from the profile are shown (i.e. within 1.5 times the electrode spacing). Finally, Fig. 9 presents pseudo 3-D models of σw and k, obtained by combining all the big-scale sections, and the 3-D geological model. The σw and k models are created by inverse distance interpolations of the 2-D sections to layers and then stacking these layers to construct a 3-D volume. In the σw 3-D model (Fig. 9a), a clear anomaly can be seen localized in the northern part of the survey, whereas in the k 3-D model (Fig. 9b) the low-permeable layer (around 10 − 12.5 to 10 − 14 m2) is more regional and represents the sand-till and lignite/clay layers (Fig. 9c). Figure 9. View largeDownload slide (a) 3-D Geological model, (b) 3-D permeability model and (c) 3-D water electrical conductivity model. Figure 9. View largeDownload slide (a) 3-D Geological model, (b) 3-D permeability model and (c) 3-D water electrical conductivity model. 5 DISCUSSION 5.1 The new approach compared to classical applications of the IP method at contaminated sites For field investigations, risk assessment and adequate remedial design at contaminated sites information on the spatial distribution of the contamination, local geology and hydrogeological properties is required. Aiming expressly at addressing these needs, we developed a new approach for characterizing contaminated sites through TD SIP data. For this purpose, the spatial distributions of two parameters directly usable in the hydrogeological interpretation were derived, that is, the distribution of water conductivity (decoupled from the surface conductivity of the medium) and hydraulic permeability. The water conductivity was used as a proxy for contamination, even though the conductivity response due to contamination is site-/process-specific (e.g. Atekwana & Atekwana 2010) and has not been targeted in this study. On the other hand, the hydraulic permeability has been used as a lithological indicator for mapping the geology (besides its intrinsic value in the hydrogeological characterization). The approach presented in this study differs significantly from the usual applications of IP for characterizing sites impacted by contamination that are generally aimed at the lithological characterization by IP parameters (e.g. Gazoty et al.2012b; Johansson et al.2016) the identification of the source of the contamination (e.g. landfill delineation, Dahlin et al.2010; Gazoty et al.2012a; Wemegah et al.2017), the identification of mobile (pools) or residual NAPLs (e.g. Orozco et al.2012; Johansson et al.2015) or of biogeochemical processes (e.g. Orozco et al.2011; Chen et al.2012). For lithological characterization, parameters related to the magnitude of polarization, such as imaginary conductivity (σ''), phase shift (ϕ) or intrinsic chargeability (m0) are generally used in IP investigations (e.g. Slater & Lesmes 2002b; Gazoty et al.2012b). However, Weller & Slater (2012) have shown that σ'' is dependent on the fluid conductivity, hence lithological interpretation based on σ'' (or ϕ and m0) might be hindered at contaminated sites if significant variability in water conductivity is present. On the contrary, as derived from laboratory results and as shown in this study, permeability estimates of unconsolidated formations depend weakly on water conductivity and hence are a better lithological indicator. Furthermore, while well-established permeability ranges are available for characterizing lithology, it is much more difficult to find appropriate ranges of IP parameters for lithological description in the literature. The proposed approach has an underlying assumption for its applicability: the contamination should have neither IP nor DC signature, except for the effect of water conductivity. This requirement is not always fulfilled, for instance in subsurface settings contaminated with NAPLs (e.g. Orozco et al.2011; Orozco et al.2012; Chen et al.2012; Johansson et al.2015). However, the IP signature is usually significant only where the contaminants are present in concentrations close to the saturation point (e.g. above 1000 mg l−1 for BTEX in Orozco et al.2012), which is generally at or close to the contamination source. Consequently, depending on the contaminant types and concentrations, the distributions of bulk conductivity and hydraulic permeability retrieved by the proposed approach can be inaccurate close to the source area. In contaminant plumes far from the source, the concentrations are usually much lower (typically μg l−1 or few mg l−1 at tens to a few hundreds of meters from the source) and the requirement of a negligible DC/IP signature is entirely fulfilled, as it is also the case at the Grindsted stream site (where for instance the maximum measured BTEX concentration was 1.7 mg l−1). Nevertheless, such plumes may pose a risk to the environment. For instance, in Denmark the limit in groundwater for benzene (a BTEX compound) is 1μg l−1. 5.2 Prediction quality and resolution The use of IP for permeability estimation in the field is not a novelty in itself (e.g. Kemna et al.2004; Hördt et al.2007; Hördt et al.2009; Attwa & Günther 2013). However, to our knowledge the expression suggested by Weller et al. (2015) for permeability estimation (eq. 4), expressly derived for unconsolidated sediments from an extensive set of samples, was used in the field only in the cross-hole survey presented by Binley et al. (2016) and in the logging-while-drilling borehole survey presented by Fiandaca et al. (2017b). Binley et al. (2016) found the field-scale IP method to be suitable for providing estimates of hydraulic permeability in coarse-grained aquifers, but to be somewhat limited for the resolution of small-scale (e.g. lenses versus layers) contrasts in hydraulic permeability variation: low contrast in permeability resulted in relatively low structural resolution based on the distribution of IP parameters. On the contrary, Fiandaca et al. (2017b) found a very good correlation (within on average one order of magnitude) between the IP-derived permeability estimates and those derived using grain size analyses and slug-tests, with similar depth-trends and permeability contrasts. The results presented in this study, obtained only from surface TD SIP measurements, almost replicate the quality of permeability prediction shown in Fiandaca et al. (2017b), except for the lower spatial resolution of the surface imaging that resulted in underestimated predictions close to geological boundaries. Indeed, the spatial resolution of the permeability prediction naturally mimics the spatial resolution of the electrical properties from which the permeability is derived. A way for improving the prediction quality of the TD SIP inversions and to mitigate the resolution issue is to incorporate prior information in the inversion process, both in terms of parameter values and correlations lengths, as done for instance in resistivity inversion by Hermans & Irving (2017). We are currently working on a model parametrization that inverts directly for hydraulic permeability, which allows for a direct integration of the hydrological prior information. The quality of the correlation between the IP-derived hydraulic permeability and the permeability derived from slug tests/grain size analyses is influenced not only by the decrease in resolution with depth of the IP imaging, but more in general by the different support volume of the two estimates. The slug tests, and even more the grain size analyses, provide very local information compared to the surface IP-based permeability estimates. Taking into account also the uncertainty in the petrophysical relationship of eq. (4), this means that a correlation stronger than what we found (Fig. 6b) is most likely not possible. Similar reasoning is valid for the correlation between the water electrical conductivity measured in groundwater and the imaged bulk conductivity (Fig. 6a). The limits in spatial resolution affect also the use of the imaging results for supporting the geological modelling: the geophysical models are sometimes smeared images of the geological models interpreted including all the borehole and geological information. Nevertheless, the imaging of both water electrical conductivity and hydraulic permeability was able to capture the main (hydro) geological units of the site and the areas of higher contamination, and gave a significant input in the site characterization. In particular, this led to the discovery of the contamination in the deep contaminated aquifer. 5.3 Applicability of petrophysical relationships The approach proposed in this study depends on petrophysical relations and it is not applicable when these relations are not valid. In particular, the relationship for permeability estimation (eq. 4) is valid only for unconsolidated, saturated samples. Saturation is also required in the laboratory-derived empirical relation used to isolate the electrolytic bulk conduction from surface conduction (eq. 2). Empirical relations have been suggested in the literature also for the prediction of permeability on consolidated sediments (e.g. Weller et al.2015) and corrections for the dependence on fluid saturation have been proposed for both bulk conductivity (Archie 1942) and surface conductivity (e.g. Vinegar & Waxman 1984; Ulrich & Slater 2004). Nevertheless, the extension of the proposed approach to unsaturated and/or consolidated media is beyond the scope of this study, which deals with saturated unconsolidated sediments. The quality and accuracy of the σbulk estimation through the BIC modelling depends on the value of l used in eq. (2) and on the assumption that the relationship is equally valid for all the investigated sediments. The good correlation found between the imaged bulk conductivity and the water conductivity measured in groundwater in this study (Fig. 6a) is not a proof of the validity of the relationship of eq. (2), but decoupling σbulk from σ'surf led to inversion models more representative of the geological understanding of the site and to reduction of the equivalence problem in the inversions. In fact, the use of the BIC model practically imposes a geometrical constraint between the imaginary conductivity ($$\sigma _{\rm{max}}^{''}$$) section and the total DC conductivity (σ0) section so that a chargeable layer is also conductive. This feature is desirable as long as the petrophysical relation between σ'surf and σ''surf is applicable, because the relation is obeyed by the inversion models by construction, while it might not be fulfilled by inversions carried out for instance with the classic Cole–Cole or the MIC modelling. 5.4 Final remarks The approach presented in this study can contribute to a cost-effective characterization of contaminated sites, reducing the number of drillings needed for the plume delineation and the definition of local geology and hydrogeological properties (also by guiding the drilling campaigns). The learnings on plume, hydraulic permeability field and geology are needed for modelling of contaminant transport, and hence for risk assessment and/or adequate remedial design. For instance, the geological model derived in this study (Fig. 9a), as well as the delineation of the contamination plume (Fig. 9c), will in later studies be used as a framework for building a numerical flow model simulation of contaminant transport towards the stream. 6 CONCLUSION We developed a new approach for characterizing contaminated sites through the imaging of water conductivity σw, hydraulic permeability k, and lithology by means of TD SIP data. We tested the approach at a contaminated stream site using 16 TD SIP profiles and an extensive set of complementary hydrologic and geologic information. For modelling the complex conductivity in the data inversion, a reparametrized Cole–Cole model was developed, namely the bulk and imaginary conductivity (BIC) model, defined in terms of the bulk conductivity σbulk, the MIC $$\sigma _{\rm {max}}^{''}$$, the relaxation time τσ, and the frequency exponent C. In the BIC model the bulk conductivity σbulk is decoupled from the real surface conductivity $$\sigma _{\rm {surf}}^{'}$$ through an empirical, laboratory-derived, petrophysical relation. The subsurface permeability distribution was estimated from the BIC inversion parameters σbulk and σ''max, using a laboratory-derived empirical equation valid for unconsolidated (saturated) sediments without any calibration. The IP-derived permeability values were found to be in good agreement with the estimates obtained using slug tests and grain size analyses: most of the estimates IP-derived k values were on average within one order of magnitude. However, for a few slug tests, made close to geological boundaries, the IP-derived k values were underestimated. This is where the smooth inversion failed to produce sharp boundaries. A positive correlation was observed between water conductivity σw measured in groundwater and the imaged bulk conductivity σbulk, and an average formation factor was estimated to F = 5.1 for converting the imaged values of bulk conductivity into water conductivity. The imaging of permeability and water conductivity allowed for a better discrimination of the clay/lignite lithology from the water conductivity, also due to the decrease in model equivalence of the inversion results, and the geophysical models were actively used for supporting the geological modelling. Furthermore, the direct comparison of the spatial distribution of the imaged water conductivity and the values measured in groundwater, in conjunction with the permeability imaging and the geological interpretation, allowed for the discovery of an unknown increase of σw in the deeper (confined) aquifer, in the northern part of the survey area. Water samples from confirmation wells drilled after the survey showed elevated concentration of inorganic compounds, which are linked to the elevated water conductivity in the confined aquifer. High concentrations of xenobiotic organic contaminants such as benzene, metabolites of chlorinated solvents, sulfonamides, and barbiturates were also observed in the new boreholes. These new findings pave the way for a detailed and inexpensive mapping of bulk/water conductivity, permeability, and lithology at contaminated sites using surface TD SIP measurements, and for cost-effective risk assessment and remedial design. Acknowledgements Support was provided by the research project GEOCON, advancing geological, geophysical and contaminant monitoring technologies for contaminated site investigation (contract 1305-00004B). The funding for GEOCON is provided by Innovation Fund Denmark. REFERENCES Aisopou A., Bjerg P.L., Sonne A.T., Balbarini N., Rosenberg L., Binning P.J., 2015. Dilution and volatilization of groundwater contaminant discharges in streams, J. Contaminant Hydrol. , 172, 71– 83. https://doi.org/10.1016/j.jconhyd.2014.11.004 Google Scholar CrossRef Search ADS   Archie G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. AIME , 146, 54– 62. https://doi.org/10.2118/942054-G Google Scholar CrossRef Search ADS   Atekwana E.A., Atekwana E.A., 2010. Geophysical signatures of microbial activity at hydrocarbon contaminated sites: a review, Surv. Geophys. , 31, 247– 283. https://doi.org/10.1007/s10712-009-9089-8 Google Scholar CrossRef Search ADS   Atekwana E.A., Atekwana E., Legall F.D., Krishnamurthy R.V., 2005. Biodegradation and mineral weathering controls on bulk electrical conductivity in a shallow hydrocarbon contaminated aquifer, J. Contaminant Hydrol. , 80, 149– 167. https://doi.org/10.1016/j.jconhyd.2005.06.009 Google Scholar CrossRef Search ADS   Attwa M., Günther T., 2013. Spectral induced polarization measurements for predicting the hydraulic conductivity in sandy aquifers, Hydrol. Earth Syst. Sci. , 17, 4079– 4094. https://doi.org/10.5194/hess-17-4079-2013 Google Scholar CrossRef Search ADS   Auken E. et al.  , 2015. An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data, Explor. Geophys. , 2015, 223– 235. https://doi.org/10.1071/EG13097 Google Scholar CrossRef Search ADS   Balbarini N., Boon W.M., Nicolajsen E., Nordbotten J.M., Bjerg P.L., Binning P.J., 2017. A 3-D numerical model of the influence of meanders on groundwater discharge to a gaining stream in an unconfined sandy aquifer, J. Hydrol. , 552, 168– 181. https://doi.org/10.1016/j.jhydrol.2017.06.042 Google Scholar CrossRef Search ADS   Barcelona M.J., 1994. Site characterization for subsurface remediation, Ground Water , 32, CONF-9410209. Basu N.B., Rao P., Poyer I.C., Annable M., Hatfield K., 2006. Flux-based assessment at a manufacturing site contaminated with trichloroethylene, J. Contaminant Hydrol. , 86, 105– 127. https://doi.org/10.1016/j.jconhyd.2006.02.011 Google Scholar CrossRef Search ADS   Benson R.C., Yuhr L.B., 2016a. Hydrologic characterization and measurements, in Site Characterization in Karst and Pseudokarst Terraines: Practical Strategies and Technology for Practicing Engineers, Hydrologists and Geologists , pp. 275– 293, eds Benson R.C., Yuhr L.B., Springer Netherlands. Google Scholar CrossRef Search ADS   Benson R.C., Yuhr L.B., 2016b. What is site characterization, in Site Characterization in Karst and Pseudokarst Terraines: Practical Strategies and Technology for Practicing Engineers, Hydrologists and Geologists , pp. 99– 106, eds Benson R.C., Yuhr L.B. Springer Netherlands. Google Scholar CrossRef Search ADS   Binley A., 2015. Tools and Techniques: DC Electrical Methods, in Treatise on Geophysics , 2nd edn, Vol. 11, pp. 233– 259, ed. Schubert G., Elsevier. Google Scholar CrossRef Search ADS   Binley A., Slater L.D., Fukes M., Cassiani G., 2005. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstone, Water Resour. Res. , 41, 1– 13. https://doi.org/10.1029/2005WR004202 Google Scholar CrossRef Search ADS   Binley A., Keery J., Slater L., Barrash W., Cardiff M., 2016. The hydrogeologic information in cross-borehole complex conductivity data from an unconsolidated conglomeratic sedimentary aquifer, Geophysics , 81, E409– E421. https://doi.org/10.1190/geo2015-0608.1 Google Scholar CrossRef Search ADS   Bjerg P.L., Tuxen N., Reitzel L.A., Albrechtsen H.-J., Kjeldsen P., 2011. Natural attenuation processes in landfill leachate plumes at three Danish sites, Ground Water , 49, 688– 705. https://doi.org/10.1111/j.1745-6584.2009.00613.x Google Scholar CrossRef Search ADS PubMed  Börner F.D., Schopper J.R., Weller A., 1996. Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements1, Geophys. Prospect. , 44, 583– 601. https://doi.org/10.1111/j.1365-2478.1996.tb00167.x Google Scholar CrossRef Search ADS   Bouwer H., Rice R., 1976. A slug test for determining hydraulic conductivity of unconfined aquifers with completely or partially penetrating wells, Water Resour. Res. , 12, 423– 428. https://doi.org/10.1029/WR012i003p00423 Google Scholar CrossRef Search ADS   Cameron R.E., 1992. Aguide for site and soil description in hazardouswaste site characterization, in Superfund Risk Assessment in Soil Contamination Studies , pp. 3– 17, ed. Hoddinott K.B., ASTM International. Chambers J.E., Kuras O., Meldrum P.I., Ogilvy R.D., Hollands J., 2006. Electrical resistivity tomography applied to geologic, hydrogeologic, and engineering investigations at a former waste-disposal site, Geophysics , 71, B231– B239. https://doi.org/10.1190/1.2360184 Google Scholar CrossRef Search ADS   Chen J., Hubbard S.S., Williams K.H., Orozco A.F., Kemna A., 2012. Estimating the spatiotemporal distribution of geochemical parameters associated with biostimulation using spectral induced polarization data and hierarchical Bayesian models, Water Resour. Res. , 48, 1– 25. https://doi.org/10.1029/2011WR010992 Google Scholar CrossRef Search ADS   Cole K.S., Cole R.H., 1941. Dispersion and absorption in dielectrics I. Alternating current characteristics, J. Chem. Phys. , 9, 341– 351. https://doi.org/10.1063/1.1750906 Google Scholar CrossRef Search ADS   Dahlin T., 2001. The development of DC resistivity imaging techniques, Comput. Geosci. , 27, 1019– 1029. https://doi.org/10.1016/S0098-3004(00)00160-6 Google Scholar CrossRef Search ADS   Dahlin T., Zhou B., 2006. Multiple-gradient array measurements for multichannel 2D resistivity imaging, Near Surf. Geophys. , 4, 113– 123. Dahlin T., Rosqvist H., Leroux V., 2010. Resistivity—IP mapping for landfill applications, First Break , 28, 101– 105. Deryck S.M., Redman J.D., Annan A.P., 1993. Geophysical monitoring of a controlled kerosene spill, in Symposium on the Application of Geophysics to Engineering and Environmental Problems 1993 , pp. 5–19, doi:10.4133/1.2922041. Devlin J.F., 2015. HydrogeoSieveXL: an Excel-based tool to estimate hydraulic conductivity from grain-size analysis, Hydrogeol. J. , 23, 837– 844. https://doi.org/10.1007/s10040-015-1255-0 Google Scholar CrossRef Search ADS   Ellis P.A., Rivett M.O., 2007. Assessing the impact of VOC-contaminated groundwater on surface water at the city scale, J. Contaminant Hydrol. , 91, 107– 127. https://doi.org/10.1016/j.jconhyd.2006.08.015 Google Scholar CrossRef Search ADS   EPA, 1991. Seminar Publication: Site Characterization for Subsurface Remediation , U.S. Environmental Protection Agency. Fiandaca G., Auken E., Gazoty A., Christiansen A.V., 2012. Time-domain-induced polarization: Full-decay forward modeling and 1D laterally constrained inversion of Cole-Cole parameters, Geophysics , 77, E213– E225. https://doi.org/10.1190/geo2011-0217.1 Google Scholar CrossRef Search ADS   Fiandaca G., Ramm J., Binley A., Gazoty A., Christiansen A.V., Auken E., 2013. Resolving spectral information from time domain induced polarization data through 2-D inversion, Geophys. J. Int. , 192, 631– 646. https://doi.org/10.1093/gji/ggs060 Google Scholar CrossRef Search ADS   Fiandaca G., Christiansen A.V., Auken E., 2015. Depth of Investigation for Multi-Parameters Inversions , pp. 666– 670, European Association of Geoscientists & Engineers Publications B.V. (EAGE). Fiandaca G., Madsen L.M., Maurya P.K., 2017a. Re-parameterization of the Cole-Cole Model for Improved Spectral Inversion of Induced Polarization Data, in 23rd European Meeting of Environmental and Engineering Geophysics , doi:10.3997/2214-4609.201702026. Fiandaca G., Maurya P.K., Balbarini N., Hördt A., Christiansen A.V., Foged N., Bjerg P.L., Auken E., 2017b. Hydraulic permeability estimation directly from logging-while-drilling induced polarization data, in AGU-SEG Hydrogeophysics Workshop , Stanford, CA. Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., 2012. Mapping of landfills using time-domain spectral induced polarization data: The Eskelund case study, Near Surf. Geophys. , 10, 575– 586. https://doi.org/10.3997/1873-0604.2012046 Google Scholar CrossRef Search ADS   Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., Pedersen J.K., 2012b. Application of time domain induced polarization to the mapping of lithotypes in a landfill site, Hydrol. Earth Syst. Sci. , 16, 1793– 1804. Google Scholar CrossRef Search ADS   Gazoty A., Fiandaca G., Pedersen J., Auken E., Christiansen A.V., 2013. Data repeatability and acquisition techniques for time-domain spectral induced polarization, Near Surf. Geophys. , 11, 391– 406. https://doi.org/10.3997/1873-0604.2013013 Google Scholar CrossRef Search ADS   Günther T., Martin T., 2016. Spectral two-dimensional inversion of frequency-domain induced polarization data from a mining slag heap, J. Appl. Geophys. , 2016, 10. Hermans T., Irving J., 2017. Facies discrimination with ERT using a probabilistic methodology: effect of sensitivity and regularization, Near Surf. Geophys. , 15( 1), 13– 25. Heron G., Bjerg P.L., Gravesen P., Ludvigsen L., Christensen T.H., 1998. Geology and sediment geochemistry of a landfill leachate contaminated aquifer (Grindsted, Denmark), J. Contaminant Hydrol. , 29, 301– 317. https://doi.org/10.1016/S0169-7722(97)00028-4 Google Scholar CrossRef Search ADS   Hinsby K., Bjerg P.L., Andersen L.J., Skov B., Clausen E.V., 1992. A mini slug test method for determination of a local hydraulic conductivity of an unconfined sandy aquifer, J. Hydrol. , 136, 87– 106. https://doi.org/10.1016/0022-1694(92)90006-H Google Scholar CrossRef Search ADS   Hördt A., Hanstein T., Hönig M., Neubauer F.M., 2006. Efficient spectral IP-modelling in the time domain, J. Appl. Geophys. , 59, 152– 161. https://doi.org/10.1016/j.jappgeo.2005.09.003 Google Scholar CrossRef Search ADS   Hördt A., Blaschek R., Kemna A., Zisser N., 2007. Hydraulic conductivity estimation from induced polarisation data at the field scale — the Krauthausen case history, J. Appl. Geophys. , 62, 33– 46. https://doi.org/10.1016/j.jappgeo.2006.08.001 Google Scholar CrossRef Search ADS   Hördt A., Druiventak A., Blaschek R., Binot F., Kemna A., Kreye P., Zisser N., 2009. Case histories of hydraulic conductivity estimation with induced polarization at the field scale, Near Surf. Geophysics. , 7, 529– 545. https://doi.org/10.3997/1873-0604.2009035 Google Scholar CrossRef Search ADS   Houmark-Nielsen M., 2007. Extent and age of Middle and Late Pleistocene glaciations and periglacial episodes in southern Jylland, Denmark, Bull. Geol. Soc. Denmark , 55, 9– 35. Houmark-Nielsen M., 2011. Pleistocene glaciations in denmark: a closer look at chronology, ice dynamics and landforms, in Developments in Quaternary Sciences , pp. 47– 58, eds Ehlers J., Gibbard P.L., Hughes P.D., Elsevier. Google Scholar CrossRef Search ADS   Hvorslev M., 1951. Time lag and soil permeability in groundwater observations, US Army Corps Eng. Waterways Exp. Sta. Bull. , 36. Hönig M., Tezkan B., 2007. 1D and 2D Cole-Cole-inversion of time-domain induced-polarization data, Geophys. Prospect. , 55, 117– 133. https://doi.org/10.1111/j.1365-2478.2006.00570.x Google Scholar CrossRef Search ADS   Johansson S., Fiandaca G., Dahlin T., 2015. Influence of non-aqueous phase liquid configuration on induced polarization parameters: Conceptual models applied to a time-domain field case study, J. Appl. Geophys. , 123, 295– 309. https://doi.org/10.1016/j.jappgeo.2015.08.010 Google Scholar CrossRef Search ADS   Johansson S., Sparrenbom C., Fiandaca G., Lindskog A., Olsson P.-I., Dahlin T., Rosqvist H., 2017. Investigations of a Cretaceous limestone with spectral induced polarization and scanning electron microscopy, Geophys. J. Int. , 208, 954– 972. https://doi.org/10.1093/gji/ggw432 Google Scholar CrossRef Search ADS   Jørgensen F., Møller R.R., Nebel L., Jensen N., Christiansen A.V., Sandersen P., 2013. A method for cognitive 3D geological voxel modelling of AEM data, Bull. Eng. Geol. Environ. , 72, 421– 432. https://doi.org/10.1007/s10064-013-0487-2 Google Scholar CrossRef Search ADS   Junejo S.A., Zhou Q.Y., Talpur M.A., Debao L., Shaikh S.A., 2015. Imaging of contaminant plumes using ERT in Qinhuai river water and its bed caused by urban effluents at Nanjing, China, Environ. Earth Sci. , 74, 7431– 7440. https://doi.org/10.1007/s12665-015-4728-5 Google Scholar CrossRef Search ADS   Kemna A., Binley A., Slater L., 2004. Crosshole IP imaging for engineering and environmental applications, Geophysics , 69, 97– 107. https://doi.org/10.1190/1.1649379 Google Scholar CrossRef Search ADS   Kemna A. et al.  , 2012. An overview of the spectral induced polarization method for near-surface applications, Near Surf. Geophys. , 10, 453– 468. https://doi.org/10.3997/1873-0604.2012027 Kemna A., Huisman J.A., Zimmermann E., Martin R., Zhao Y., Treichel A., Flores Orozco A., Fechner T., 2014. Broadband electrical impedance tomography for subsurface characterization using improved corrections of electromagnetic coupling and spectral regularization, in Tomography of the Earth's Crust: From Geophysical Sounding to Real-Time Monitoring: GEOTECHNOLOGIEN Science Report No. 21 , pp. 1– 20, eds Weber M., Münch U., Springer International Publishing. Google Scholar CrossRef Search ADS   Lesmes D., Frye K.M., 2001. Influence of pore fluid chemistry on the complex conductivity and induced polarization responses of Berea sandstone, J. geophys. Res. , 160, 4079– 4090. https://doi.org/10.1029/2000JB900392 Google Scholar CrossRef Search ADS   Madsen L.M., Kirkegaard C., Fiandaca G., Christiansen A.V., Auken E., 2016. An analysis of Cole–Cole parameters for IP data using Markov chain Monte Carlo, in IP2016—4th International Workshop on Induced Polarization , Aarhus, Denmark. Madsen L.M., Fiandaca G., Christiansen A.V., Auken E., 2018. Resolution of well-known resistivity equivalences by inclusion of time-domain induced polarization data, Geophysics , 83, E47– E54. Google Scholar CrossRef Search ADS   Maurya P.K., Rønde V.K., Fiandaca G., Balbarini N., Auken E., Bjerg P.L., Christiansen A.V., 2017. Detailed landfill leachate plume mapping using 2D and 3D electrical resistivity tomography—with correlation to ionic strength measured in screens, J. Appl. Geophys. , 138, 1– 8. https://doi.org/10.1016/j.jappgeo.2017.01.019 Google Scholar CrossRef Search ADS   Olsson P.I., Dahlin T., Fiandaca G., Auken E., 2015. Measuring time-domain spectral induced polarization in the on-time:decreasing acquisition time and increasing signal-to-noise ratio, J. Appl. Geophys. , 123, 316– 321. Google Scholar CrossRef Search ADS   Olsson P.-I., Fiandaca G., Larsen J.J., Dahlin T., Auken E., 2016. Doubling the spectrum of time-domain induced polarization by harmonic de-noising, drift correction, spike removal, tapered gating and data uncertainty estimation, Geophys. J. Int. , 207, 774– 784. https://doi.org/10.1093/gji/ggw260 Google Scholar CrossRef Search ADS   Orozco A.F., Williams K.H., Long P.E., Hubbard S.S., Kemna A., 2011. Using complex resistivity imaging to infer biogeochemical processes associated with bioremediation of an uranium-contaminated aquifer, J. geophys. Res. , 116, G03001, doi:10.1029/2010JG001591. Orozco A.F., Kemna A., Oberdörster C., Zschornack L., Leven C., Dietrich P., Weiss H., 2012. Delineation of subsurface hydrocarbon contamination at a former hydrogenation plant using spectral induced polarization imaging, J. Contaminant Hydrol. , 136–137, 131– 144. https://doi.org/10.1016/j.jconhyd.2012.06.001 Google Scholar CrossRef Search ADS   Panagos P., Van Liedekerke M., Yigini Y., Montanarella L., 2013. Contaminated sites in Europe: review of the current situation based on data collected through a European network, J. Environ. Public Health , 2013, 1– 11. https://doi.org/10.1155/2013/158764 Google Scholar CrossRef Search ADS   Pelton W.H., Ward S.H., Hallof P.G., Sill W.R., Nelson P.H., 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP, Geophysics , 43, 588– 609. https://doi.org/10.1190/1.1440839 Google Scholar CrossRef Search ADS   Rasmussen E.S., Dybkjær K., Piasecki S., 2010. Lithostratigraphy of the upper Oligocene - Miocene succession of Denmark, Geol. Surv. Denmark Greenland Bull. , 22, 1– 92. Rasmussen J.J., McKnight U.S., Sonne A.T., Wiberg-Larsen P., Bjerg P.L., 2016. Legacy of a chemical factory site: contaminated groundwater impacts stream macroinvertebrates, Arch. Environ. Contam. Toxicol. , 70, 219– 230. https://doi.org/10.1007/s00244-015-0211-2 Google Scholar CrossRef Search ADS PubMed  Revil A., 2012. Spectral induced polarization of shaly sands: Influence of the electrical double layer, Water Resour. Res. , 48, W02517, doi:10.1029/2011WR011260. https://doi.org/10.1029/2011WR011260 Revil A., Florsch N., 2010. Determination of permeability from spectral induced polarization in granular media, Geophys. J. Int. , 181, 1480– 1498. Rønde V., McKnight U.S., Sonne A.T., Balbarini N., Devlin J.F., Bjerg P.L., 2017. Contaminant mass discharge to streams: Comparing direct groundwater velocity measurements and multi-level groundwater sampling with an in-stream approach, J. Contaminant Hydrol. , 206, 43– 54. https://doi.org/10.1016/j.jconhyd.2017.09.010 Google Scholar CrossRef Search ADS   Royse K.R., 2010. Combining numerical and cognitive 3D modelling approaches in order to determine the structure of the Chalk in the London Basin, Comput. Geosci. , 36, 500– 511. https://doi.org/10.1016/j.cageo.2009.10.001 Google Scholar CrossRef Search ADS   Sims R.C., 1990. Soil remediation techniques at uncontrolled hazardous waste sites, J. Air Waste Manage. Assoc. , 40, 704– 732. https://doi.org/10.1080/10473289.1990.10466716 Google Scholar CrossRef Search ADS PubMed  Slater L., 2007. Near surface electrical characterization of hydraulic conductivity: from petrophysical properties to aquifer geometries-A Review, Surv. Geophys. , 28, 169– 197. https://doi.org/10.1007/s10712-007-9022-y Google Scholar CrossRef Search ADS   Slater L., Lesmes D.P., 2002a. Electrical-hydraulic relationships observed for unconsolidated sediments, Water Resour. Res. , 38, 31-1–31-13. https://doi.org/10.1029/2001WR001075 Slater L.D., Lesmes D., 2002b. IP interpretation in environmental investigations, Geophysics , 67, 77– 88. https://doi.org/10.1190/1.1451353 Google Scholar CrossRef Search ADS   Smith S.H., 1962. Temperature correction in conductivity measurements, Limnol. Oceanogr. , 7, 330– 334. https://doi.org/10.4319/lo.1962.7.3.0330 Google Scholar CrossRef Search ADS   Sonne A.T., McKnight U.S., Rønde V.K., Bjerg P.L., 2017. Assessing the chemical contamination dynamics in a mixed land use stream system, Water Res ., 125, 141– 151. https://doi.org/10.1016/j.watres.2017.08.031 Google Scholar CrossRef Search ADS PubMed  Sparrenbom C.J., Åkesson S., Johansson S., Hagerberg D., Dahlin T., 2017. Investigation of chlorinated solvent pollution with resistivity and induced polarization, Sci. Total Environ. , 575, 767– 778. https://doi.org/10.1016/j.scitotenv.2016.09.117 Google Scholar CrossRef Search ADS PubMed  Springer R., Gelhar L., 1991. Characterization of large-scale aquifer heterogeneity in glacial outwash by analysis of slug tests with oscillatory response, Cape Cod, Massachusetts, US Geol. Surv. Water Res. Invest. Rep. , 91, 36– 40. Switzer A.D., Pile J., 2015. Grain size analysis, in Handbook of Sea-Level Research , eds Shennan I., Long A. J., Horton B. P., John Wiley & Sons, p. 331. Google Scholar CrossRef Search ADS   Sørensen K.I., Larsen F., 1999. Ellog auger drilling: three-in-one method for hydrogeological data collection, Ground Water Monitor Remediation , 19, 97– 101. https://doi.org/10.1111/j.1745-6592.1999.tb00245.x Google Scholar CrossRef Search ADS   Tarasov A., Titov K., 2013. On the use of the Cole–Cole equations in spectral induced polarization, Geophys. J. Int. , 195, 352– 356. https://doi.org/10.1093/gji/ggt251 Google Scholar CrossRef Search ADS   Ulrich C., Slater L.D., 2004. Induced polarization measurements on unsaturated, unconsolidated sands, Geophysics , 69, 762– 771. https://doi.org/10.1190/1.1759462 Google Scholar CrossRef Search ADS   Vinegar H.J., Waxman M.H., 1984. Induced polarization of shaly sands, Geophysics , 49, 1267– 1287. https://doi.org/10.1190/1.1441755 Google Scholar CrossRef Search ADS   Weller A., Slater L., 2012. Salinity dependence of complex conductivity of unconsolidated and consolidated materials: Comparisons with electrical double layer models, Geophysics , 77, D185– D198. https://doi.org/10.1190/geo2012-0030.1 Google Scholar CrossRef Search ADS   Weller A., Breede K., Slater L., Nordsiek S., 2011. Effect of changing water salinity on complex conductivity spectra of sandstones, Geophysics , 76, F315– F327. https://doi.org/10.1190/geo2011-0072.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Nordsiek S., 2013. On the relationship between induced polarization and surface conductivity: implications for petrophysical interpretation of electrical measurements, Geophysics , 78, D315– D325. https://doi.org/10.1190/geo2013-0076.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Binley A., Nordsiek S., Xu S., 2015. Permeability prediction based on induced polarization: insights from measurements on sandstone and unconsolidated samples spanning a wide permeability range, Geophysics , 80, D161– D173. https://doi.org/10.1190/geo2014-0368.1 Google Scholar CrossRef Search ADS   Wemegah D.D., Fiandaca G., Auken E., Menyeh A., Danuor S.K., 2017. Spectral time-domain induced polarisation and magnetic surveying—an efficient tool for characterisation of solid waste deposits in developing countries, Near Surf. Geophys. , 15, 75– 84. Yang C.H., Yu C.Y., Su S.W., 2007. High resistivities associated with a newly formed LNAPL plume imaged by geoelectric techniques—a case study, J. Chin. Inst. Eng. , 30, 53– 62. https://doi.org/10.1080/02533839.2007.9671230 Google Scholar CrossRef Search ADS   Zisser N., Kemna A., Nover G., 2010. Dependence of spectral-induced polarization response of sandstone on temperature and its relevance to permeability estimation, J. geophys. Res. , 115, doi:10.1029/2010JB007526. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations