Spectral induced polarization of partially saturated clay-rocks: a mechanistic approach

We have developed a mechanistic model to interpret spectral induced polarization data of partially saturated clay-rocks. This model accounts for the polarization of the grains through an electrical double layer model with a polarization model of the inner part of the electrical double layer called the Stern layer. The polarization model accounts also for the Maxwell–Wagner polarization at frequencies higher than 100 Hz. The Maxwell–Wagner polarization is modelled by using a conductivity model modiﬁed to account for the presence of a non-wetting immiscible phase like air in the pore space. The resulting model is consistent with the ﬁrst and second Archie’s laws in the case where surface conductivity can be neglected. The volumetric charge density of the diffuse layer at saturation is divided by the saturation of the water phase to account for the partial water saturation of the porous material. The model comprises seven fundamental parameters: the formation factor, the second Archie’s exponent, a critical water saturation level, the mean electrical potential of the pore space at saturation, the density of the counterions in the Stern layer, and at least two parameters describing the grain size distribution. Most of these parameters can be derived independently using alternative measurements and electrochemical models. Measurements were performed in the frequency range 10 mHz– 45 kHz using ﬁve samples from the Callovo–Oxfordian formation in the eastern part of the Paris Basin, France. The model agrees fairly well with the experimental data at saturation and for partially saturated clay-rocks down to 1 Hz. Most of the seven physical parameters entering the model were independently evaluated.


I N T RO D U C T I O N
Spectral induced polarization reports the modulus of the conductivity and the phase between the current and the voltage in the frequency domain (typically from few millihertz to few tens of kilohertz, sometimes few Megahertz). The modulus and the phase can also be recasted into an equivalent complex conductivity or complex resistivity. This method has been developed first to locate ore bodies (Marshall & Madden 1959;Van Voorhis et al. 1973) and more recently in environmental sciences (Kemna et al. 2000(Kemna et al. , 2004Slater & Binley 2006;Hördt et al. 2007;Blaschek et al. 2008). In environmental sciences, spectral induced polarization has been used to map contaminants in the ground (Vanhala 1997b), to determine grain size distribution (Lesmes & Morgan 2001;Leroy et al. 2008) and permeability (Börner et al. 1996;de Lima & Niwas 2000;Binley et al. 2005;Hördt et al. 2007). In parallel to these new types of investigations, the improvement of high accuracy acquisition systems able to measure the phase with a precision of 1-0.1 mrad (Zimmermann et al. 2008) has allowed to investigate rock samples and soils characterized by very low quadrature conductivity including glass beads (Leroy et al. 2008), clean (clay-free), sandstones (Binley et al. 2005) and clay materials ).
During the last few years, several works have been done in using spectral induced polarization to study the properties of clay-rocks (Kruschwitz & Yaramanci 2004;Cosenza et al. 2008;Tabbagh et al. 2009). These works were motivated by the use of the spectral induced polarization method to study nuclear waste disposals (Kruschwitz & Yaramanci 2004;Pusch 2006) and generally speaking to study the role of clay barriers for the confinement of different types of hazardous materials. Indeed, because of their high-specific surface area and very low permeabilities (<10 −19 m 2 ) (Gaucher et al. 2004), clay-rocks are good candidates for host rocks (Homand Spectral induced polarization of clay-rocks 211 et al. 2004). Various national agencies (like ANDRA in France) have been commissioned to study the transport properties, speciation and thermohydro-mechanical behaviour of clay-rocks. These works will allow to assess the long-term safety of future potential radioactive waste disposals in these materials.
In order to study these clay-rocks formations, underground research laboratories have been built like Bure in France and Mol in Belgium or test sites have been selected, like Tournemire in France, and Mont Terry in Switzerland. Gallery excavations in these laboratories have generated perturbations of the rock: the Excavated Damaged Zone (EDZ) along the walls of the galleries (Tsang et al. 2005). To improve the knowledge of these zones, geophysical methods can provide very useful information especially spectral induced polarization as shown by Kruschwitz & Yaramanci (2004). Additional investigations were performed in the laboratory, which confirm this possibility Ghorbani et al. 2009). However, all these authors used semi-empirical approaches based on the Cole-Cole model (Cole & Cole 1941;Ghorbani et al. 2007;Chen et al. 2008) or similar semi-empirical models like the Davidson-Cole model (Davidson & Cole 1950) and the generalized Cole-Cole model (Vanhala 1997a;Ghorbani et al. 2009).  have proposed recently a mechanistic model to quantify spectral induced polarization in saturated clay-rich materials. However, a quantitative understanding of spectral induced polarization data in unsaturated clay-rocks is still missing.
We propose below a mechanistic model of spectral induced polarization of unsaturated clays-rocks. This model is based on an extension of the theory developed recently by  for clay-rich materials. This model extends also some previous modelling efforts to understand cross-coupling properties of clay materials Revil & Linde 2006). We also compare the prediction of this new model to new spectral induced polarization measurements performed on five core samples from the Callovo-Oxfordian clay-rock formation. This formation is presently studied by ANDRA, the French Nuclear Waste Agency, as a potential host for nuclear wastes.

D E S C R I P T I O N O F T H E C A L L OV O -O X F O R D I A N C L AY-RO C K
ANDRA has developed an underground research facility, the test site of Bure, to study the possibility of long-term underground disposal of nuclear wastes. The site of Bure (Meuse/Haute-Marne, France) is located in the eastern part of the Paris Basin in a thick clay-stone formation of Callovian to Oxfordian age ( Fig. 1). At Bure, the Callovo-Oxfodian formation has a thickness of approximately 130-160 m. The Callovo-Oxfordian formation can be subdivided in four main lithofacies units. At Borehole EST205 (see location in Fig. 1), we find first the C2d unit for depths comprised between approximately 417 and 437 m. This unit presents some smooth spatial variations in the mineralogy and petrophysical properties. The C2c and C2b units are located between approximately 437 and 456 m. These units are relatively homogeneous. The bottom unit C2b of the Callovo-Oxfordian formation unit is itself divided into two subunits: C2b2 (approximately 456-473 m) and C2b1 (approximately 473-508 m). The subunit C2b1 corresponds to the unit where the proportion of clay minerals is the highest through the Callovo-Oxfordian formation and is located at depths comprised between approximately 486 and 489 m. The underground laboratory of ANDRA at Bure is located just below this level. The bottom unit   of the Callovo-Oxfordian formation corresponds to the C2a unit (<508 m). The mean gravimetric water contents and porosities of the different units are summarized in Table 1 (ANDRA 2005). The Callovo-Oxfordian clay-rock can be conceptualized as a clay matrix with some imbedded grains of silica and carbonates (calcite and a minor fraction of dolomite). (Fig. 2). The clay matrix represents 20-50 per cent of the rock volume in the total formation (up to 45-50 per cent in the subunit C2b1). This clay fraction is mainly composed of illite and interstratified illite-smectite clays. There is also a small amount of kaolinite in the lower part of the formation. Smectite based clay-minerals form aggregates called tactoids. The size of the tactoid is typically 10 μm (ANDRA 2005). The size of the silica and carbonate grains is usually larger than 4 μm with a mean around 10-20 μm. The volume of the silica and carbonate grains represent between 20 and 40 per cent of the rock assemblage. In addition, there are small amounts of potassic feldspar, plagioclase and pyrite. An extensive geochemical analysis of the Callovo-Oxfordian formation has been proposed recently by Gaucher et al. (2004).

Theory
We consider a representative elementary volume of the Callovo-Oxfordian clay-rock of porosity φ. The pore volume is  (Gaucher et al. 2006) (only charged species) a Obtained using PHREEQC2 (see Leroy et al. 2007). The software PHREEQC2 is described in Parkhurst & Appelo (1999).
filled by two continuous and immiscible fluid phases, a wetting fluid (subscript w) and a non-wetting fluid. In our analysis, the wetting fluid is water and the non-wetting fluid is air. We denote s w ∈ [0; 1] the saturation of the water phase. This saturation is related to the volumetric water content θ by s w = θ/φ. In soil science, many authors define also a gravimetric water content w defined as the ratio between the mass of the water phase m w relative to the mass of the solid phase m s of the core sample. Following Revil et al. (2007) and , we define a percolation threshold s w c for the water phase. Below this threshold, the water phase is not continuous at the scale of the representative elementary volume (see also Hunt & Ewing 2003 for the use of percolation theory in the context of this problem).
The minerals of the Callovo-Oxfordian clay-rock, especially the clay minerals, have negatively charged surface at the pH of the pore water (pH = 7.2, see Table 2; Leroy et al. 2007). To counterbalance this negative charge on the mineral surface, the pore water has an excess of cations (see Leroy et al. 2007 for a geochemical modelling of the pore water composition of the Callovo-Oxfordian clay-rock). These counterions are partitioned between the Stern layer, which is attached to the mineral surface, and the diffuse layer located inside the pore water phase (Fig. 3). In our modelling approach and in unsaturated condition, the Callovo-Oxfordian clay-rock is composed by three phases: (1) the solid phase, including the insulating minerals coated by the conductive Stern layer, (2) the water phase including the diffuse layer and locally some free electrolyte pockets in some macropores and (3) the non-wetting phase corresponding to air. Neglecting the surface charge of the air-water interface (S aw ) (see Linde et al. 2007;Revil et al. 2007), the charge balance  Figure 3. Sketch of the distribution of the ionic species in the pore space of a charged unsaturated porous medium. We note S sw the shear plane between the mineral surface and the pore water and S wa the interface between the air phase and the pore water (these two interfaces have complex shapes). The pore water is characterized by a volumetric charge densityQ V . ϕ 0 , ϕ β and ϕ d are the electrical potentials of the surface of the mineral surface, the midplane of the β-Stern layer, and the shear plane, respectively. When the mean distance between the air/water interface and the mineral/water interface decreases, the ions of the pore water phase are confined in a smaller volume and therefore the conductivity of the pore water phase increases.
condition is where Q S is the total surface charge density (in Cm −2 ) at the solidsolution interface of area S sw (including the true charge density on the mineral surface plus the surface charge density of the Stern layer), V w (in m 3 ) denotes the volume of the water phase andQ V /s w represents the excess of charge per unit of volume of water in the unsaturated porous material (Q V is the volumetric charge density at saturation of the water phase). Considering a solution with Q ionic species i, the effective charge density of the pore water is written as where q i is the charge (± z i e, where z i is the valency of the ion and e the charge of the electron, that is, ∼1.6 × 10 −19 C) andC i the ionic concentration of the ion i in the medium bulk solution (in m −3 ). The bulk pore water solution of the clay matrix is mainly composed by the diffuse layer due to the small size of the pores (Leroy et al. 2007;Jougnot et al. 2009). We can use the Donnan approach to determine the ionic concentrations in the pore space of the clay matrixC i (Revil & Linde 2006;Leroy et al. 2007;Revil 2007).
The ionic species that are present in the pore space of the clayrock are assumed to be in local thermodynamic equilibrium with an infinite reservoir of the same Q ionic species at concentrations C i .
Equilibrium implies (Revil & Linde 2006) where ϕ m represents the mean electrical potential in the pore water phase (in V), k B represents the Boltzmann constant (k B = 1.3806 × 10 −23 JK −1 ) and T is the absolute temperature (in K). A very important parameter in the study of the electrical double layer is the partition coefficient f Q Leroy et al. 2007). This partition coefficient represents the fraction of the countercharge located in the Stern layer to compensate the true surface charge density of the mineral surface. We note 0 Xi the surface site density of species i in the Stern layer and d i the surface site density of species i in the diffuse layer (both in m −2 ). The partition coefficient f Q is defined as (Leroy et al. 2007).
where z i is the charge number of species i. For unsaturated conditions, the total charge density of counterions per unit of pore volume, Q V , can be related to the excess of charge per unit of volume of solution in the medium,Q V /s w and to the partition coefficient defined above bȳ and Q V can be related in turn to the cation exchange capacity (CEC, in C kg −1 ) of the rock: where ρ s is the bulk density of the solid phase of the material. The CEC characterizes the density of charge on the mineral surface because it can be measured directly using chemical titration methods. The electrical double layer theory developed by Leroy et al. (2007) includes a specific speciation model of the different crystalline planes of a mineral. This model can be used to infer the dependence of the partition coefficient f Q with the ionic strength of the pore water, Leroy  proposed recently a model to study the spectral induced polarization of clay-rich porous materials saturated by water. This model is based on two different polarization mechanisms: (1) a polarization of the Stern layer associated with the discontinuity of the solid phase and (2) a Maxwell-Wagner polarization process. This model has been built for the case of a binary symmetric electrolyte (like NaCl or KCl). This is a major assumption of this study as discussed further below. Considering the composition of the pore water proposed by Gaucher et al. (2006) and Leroy et al. (2007) (see Table 2), we see that the sodium Na + is the dominating counterion of the pore water. Therefore using a model with a single type of counterions is a reasonable assumption. However other counterions may be present and may play some role in widening the spectrum of the relaxation times. The sodium cation is located on the Outer Helmoltz Plane of the Stern layer, which is also called the d-plane (see Fig. 3), because sodium keeps its hydration shell on the mineral surface (Tournassat et al. 2009). Therefore it keeps some mobility along the mineral surface. If sodium would form an inner sphere complex with the mineral surface, it would have no tangential mobility along the mineral surface.

D. Jougnot et al.
To extend the spectral induced polarization model developed by Leroy et al. (2007) to unsaturated conditions, and following Revil et al. (2007), we use the following change of variables to introduce the effect of the saturation of the water phase: Revil 1999), where F is the formation factor and n is the saturation exponent also called the second Archie's exponent (Archie 1942).
The polarization of the Stern layer is obtained by using the approach developed by Schwarz (1962). This approach is based on the assumption that there are no exchange of ions between the Stern and the diffuse layers. This assumption seems valid because the kinetics of ion sorbtion/desorption in the Stern layer is slow (typically several hours). We consider therefore that the movement of the ions in the Stern layer follows tangentially the pore water/mineral interface and that the characteristic length scale associated with this polarization mechanism is the size of the grains. In the model developed by Schwarz (1962), the relaxation time for the electromigration of an ion of species i, τ i (in s), is therefore related to the particle radius a (in m) and to the diffusion coefficient of the ion i in the Stern layer D i S by, This time constant corresponds to the relaxation time of the diffusion process of the counterions in the Stern layer with a characteristic distance a. In our model, the diffuse layer is homogeneized with the pore water phase using the Donnan equilibrium approach (Revil & Linde 2006). We do not consider polarization for the diffuse layer because it forms a continuous phase through the porous material. In contrast, the Stern layer polarizes because it is discontinuous.
Because the clay fraction of the Callovo-Oxfordian clay-rock is mainly composed of illite and interstratified illite-smectite clays, we will consider only the response of illite and smectite to define the spectral induced polarization of the clay-rock. These minerals are coated by a Stern layer. We assume that all counterions are Na + . If all grains had the same grain radius a 0 , the Stern layer complex conductivityσ S could be computed by .
with ω the angular frequency (in rad s −1 , ω = 2πf is the frequency in Hertz), a 0 is a grain size and τ 0 Na is the associated relaxation time associated with the backdiffusion of the counterions in the Stern layer. At very low frequencies ω << 1/τ 0 Na , there is no contribution from the Stern layer to surface conductivity because of the diffusion of the counterions. At high frequencies, diffusion does not have the time to occur and therefore conductivity increases because of the contribution of the Stern layer conductivity. The parameter ∞ S corresponds therefore to the contribution of the Stern layer at high frequencies. This contribution is defined by ).
(see also Revil & Leroy 2001) where e is the elementary charge, β S Na is the ionic mobility of the sodium in the Stern layer (in m 2 s −1 V −1 ). The mobility β S Na and the diffusion coefficient D i S are related to each other by the Nernst-Einsten relationship. The complex conductivity of the Stern layer σ * S (in Sm −1 ) is determined as ).
where ε s represents the dielectric permittivity of the clay minerals (ε 0 = 8.85 × 10 −12 Fm −1 , the dielectric permittivity of vacuum). The effect of the grain size distribution is obtained by convoluting the response obtained for a single grain with the particle size distribution (PSD) (see ). The PSD is related to the distribution of the relaxation times through eq. (8).
The conductivity of the pore water needs to be computed accounting for the presence of Q ionic species and the effect of the diffuse layer. Generalizing the result obtained by Revil & Linde (2006) for the saturated case, the conductivity of the water phase (in Sm −1 ) in unsaturated condition is defined by, where β i is the ionic mobility of the species i in the solution (in m 2 s −1 V −1 ). We see therefore that the conductivity of the water phase depends strongly on the saturation of the water phase and increases with the desaturation. A complex pore water conductivity σ * f can be defined as, with ε f = 81ε 0 the dielectric permittivity of water (at standard atmospheric pressure and at 25 • C).
We need now an expression for the conductivity of the representative elementary volume, σ * , in unsaturated conditions. This is done by extending the result obtained by de Lima & Sharma (1992) and Revil et al. (2007) This model accounts for the Maxwell-Wagner polarization. This electrical conductivity model is part of the generalized framework proposed by Revil & Linde (2006) and Revil et al. (2007) to describe transport properties of saturated and unsaturated microporous media. This framework presents the important advantage to propose the determination of the macroscopic transport properties in a porous medium from a couple of intrinsic material properties that can be transposed from one macroscopic property to the other. From these fundamental properties, it becomes therefore possible to predict a wide variety of transport properties of the porous material (e.g. Jougnot et al. 2009 used the formation factor determined from electrical conductivity measurements to determine the diffusion coefficient of ionic tracers in core samples of the Callovo-Oxfordian clay-rock).

Sensitivity analysis
The previous model of induced polarization is based on seven independent parameters: (1) the formation factor, F, (2) the saturation index, n, (3) the percolation threshold, s c w , (4) the Donnan equilibrium potential, ϕ m , (5) the surface site density of the counterions in the Stern layer, 0 Na and (6) and (7) at least two parameters describing the grain size distribution (a mean and a standard deviation for a log normal probability distribution for instance).
The formation factor of the Callovo-Oxfordian clay-rock, in unit C2b1 and C2b2, is comprised between 15 and 140 ). The saturation exponent n has been studied by    clay-rock. They obtained n = 2.0. No values exist for s c w for the Callovo-Oxfordian clay-rock.  previously considered s c w = 0.1 as a reasonable value. We will see below that a correct estimate seems to be higher than this value. Considering the pore water composition determined by Gaucher et al. (2006) (the THERMOAR model reported in Table 2), the Donnan equilibrium potential, ϕ m , is comprised between −54 and −24 mV (see Jougnot et al. 2009). From , the surface site density, 0 Na , is comprised between 0.632 nm −2 for pure smectite to 1.12 nm −2 for pure illite. The particle size distribution in consolidated clay-stones is difficult to determine experimentally because of the existence of the tactoids. As discussed above in Section 3.1, the distribution of the relaxation times, τ , is related to the grain size distribution. We consider below a Cole-Cole distribution to characterize the distribution of the relaxation times (Cole & Cole 1941;Bötcher & Bordewijk 1978): where τ 0 , is the main relaxation time (in s) and c is the Cole-Cole exponent. This distribution is symmetric for τ = τ 0 .  obtained c = 0.76 ± 0.05 for the Mancos clay-rock and c ∈ [0.60; 0.84] for illites and smectites. This distribution is convoluted with the expression obtained for the Stern layer surface conductivity for a single grain, see eq. (9). In our model, the grain size distribution follows a Cole-Cole distribution and to this distribution is associated a distribution of relaxation time that follows also a Cole-Cole distribution. This is very different than taking an empirical Cole-Cole model to describe the modulus and the phase of the resistivity.
To perform a sensitivity analysis, the values of the seven keyparameters of our spectral induced polarization model can be chosen as follow: F = 80, n = 2.0, s c w = 0.10, ϕ m = −40 mV, 0 Na = 1 nm −2 , c = 0.85 (smectite-A in ) and τ 0 = 3.80 × 10 −2 s (this corresponds to a grain diameter equal to 10 −5 m and an ionic mobility in the Stern layer of β S Na = 5.19 × 10 −8 m 2 s −1 V −1 , according to eq. (8)). We present the spectral induced polarization data of our sensitivity test by showing the mod-ulus of the complex resistivity |ρ * | (the inverse of the complex conductivity in Ohm m) and the phase ϕ (ϕ = −arctan[Im(σ * )/Re(σ * )], in rad where Im(σ * ) represents the quadrature conductivity and Re(σ * ) represents the in phase conductivity). By convention, we will always plot (−ϕ) expressed in mrad.
The first parameter we investigate is the saturation s w . Fig. 4 shows the influence of the saturation with s w ∈ {0.2; 1}. As expected, the amplitude of the complex resistivity increases with the desaturation of the porous material because of the insulating non-wetting phase (Fig. 4a). The phase also increases when the saturation decreases (Fig. 4b). This is because the presence of a non-wetting and insulating phase creates a new interface having an effect on the Maxwell-Wagner polarization. However, because the density of the counterions has also changed in the pore water, the presence of air has also an effect on the polarization of the Stern layer. Fig. 5 shows the influence of F, n and s c w upon the spectral induced polarization response of a clay-rock. All these parameters influence the response in the same way: when they increase, both the amplitude of the complex resistivity, |ρ * | and the phase, −ϕ, increase (see Fig. 5). The sensitivity of our model to the two geochemical parameters (ϕ m and 0 Na ) is shown in Fig. 6. The Donnan equilibrium potential influences the signal in the same way as F, n and s c w . The surface site density influences the model in a completely different way. When 0 Na increases, the amplitude increases at high frequencies (>1 Hz) but not at low frequency (<1 Hz) where the phase decreases (Fig. 6). This implies that the effect of the surface site densities is different on the Maxwell-Wagner polarization, which dominates at high frequencies. This influence is made through the effect of surface conductivity, which changes the apparent conductivity of the grains. At low frequencies (below 10 Hz in the present case), the response is dominated by the polarization of the Stern layer.
We investigate now the influence of the distribution of the relaxation times by analysing the sensitivity of the model to the two parameters τ 0 and c (see Fig. 7). We note that the mean relaxation time and the Cole-Cole exponent influence the high frequency part of the amplitude (>10 Hz) (Figs 7a and d). c also modifies the phase mainly at high frequency (Fig. 8b)  lower value of c corresponds to a broader distribution of grain size. Consequently, this leads to a broader phase distribution, because the phase is strongly linked to the Stern layer polarization [the Stern layer contribution to the surface conductivity depends on the grain, size, see eq. (9)].

Description of the rock samples
In the present contribution, we only study the spectral induced polarization response of Callovo-Oxfordian clay-rock sam-ples. Due to the formation depth and the diagenetic phenomena, the Callovo-Oxfordian clay-stone is a very consolidated rok (ANDRA 2005) and samples do not fall apart after extraction. The samples used for the experiments have been extracted from the Callovo-Oxfordian layer at Bure (Meuse/Haute-Marne, France) by drilling with a NaCl saturated mud. The core samples have been extracted from the C2b1 and C2b2 units. The clay-rock has been sampled as a series of ∼20 cm long cylinders at different depths. The EST433 samples have been resampled to avoid the NaCl contamination and then became parallelepiped-like. All the samples have been conditioned in double hermetic aluminium bags and confined in polystyrene during the transport to avoid the perturbations.  transportation. While opening, the samples were already partially unsaturated (s i w ∈ [0.75; 0.81]). In order to limit the perturbation in the clay-stone pore space, we have decided not to resaturate the samples prior to our experiments. For the experiments, we used five samples from two boreholes at different depths (from units C2b1 and C2b2) and with different sizes (Table 3). We resampled some of the broken samples to make them usable for our experiments. Note that the resampling has been done with dry air to avoid sample pollution and/or osmotic perturbation (which could yield to the sample disintegration).

Experimental setup
The Spectral Induced Polarization method is based on the measurements of the complex resistivity ρ * (in m) over several decades of frequency (usually from 1 mHz to 10 kHz). The complex resistivity is the inverse of the complex conductivity:ρ * = 1/σ * . What is actually measured is the complex impedance Z * (in ) between the end-faces of a cylindrical sample for instance (Pelton et al. 1983): where U is the voltage between the end-face boundaries of the sample, I the magnitude of the current which goes through the sample, and |Z * (ω)| and ϕ(ω) are the amplitude and the phase of the complex impedance, respectively. The complex resistivity ρ * is related to Z * by a geometrical factor K (in m) by: ρ * (ω) = K Z * (ω). This geometrical factor takes into account not only the position of the electrodes on the sample, but also the size and the shape of the samples. The complex resistivity ρ * can be recasted into a modulus |ρ * (ω)| and a phase ϕ(ω). In the present contribution, amplitude will always stand for the resistivity amplitude. The phase is the same for the resistivity and the impedance. Considering the geometry of samples #EST27906, #EST29296 and #EST29300, we choose the same electrode array than proposed by Cosenza et al. (2007). It is a Wenner-α type array (see Figs 8a and c for details). This setup is called the Lg-setup below. The lengths of samples #EST28144 and #EST30749 were too short to use the Lg-setup. Therefore, we tried a different approach by putting the electrodes on the side of the core samples (name the Rd-setup below). Instead of the length of the sample, we used a constant angle α = 45 • between each electrode (see Figs 8b and d for details). Due to the orientation of the boreholes in the Callovo-Oxfordian formation, the Lg-setup corresponds to electrical measurements perpendicular to the stratification while the Rd-setup provides measurements parallel to the stratification.
Both setups are composed by two current electrodes (C1 and C2) and two potential electrodes (P1 and P2). All these electrodes were developed in the medical world. We found that they allow the best contact, and thus the lowest contact resistance, among various tested electrodes. The current (or source) electrodes are used to inject an alternative current at different frequencies. We choose the round 50 mm diameter electrodes to fit the different sample sizes. They have thin carbon films to improve the contact (Dura-Stick II). The potential electrodes are electrocardiogram (ECG) electrodes: their silver-silver chloride (Ag/AgCl) composition made them nonpolarizable (3M Red Dot). Each of these electrodes correspond to a round 10 mm diameter metal disk, galvanized with silver and covered with a soft sponge imbibed of a AgCl gel. For the computation of the geometrical constant K, it is important to note that the current electrodes have a diameter of 50 mm while the potential electrodes are nearly punctual.
The non-polarizability of the potential electrodes has been validated using porous and electrically inert samples. For example, a clean (clay-free) limestone sample saturated with a saline NaCl brine (0.78 Sm −1 at 25 • C) was used to perform a test as prescribed by Vinegar & Waxman (1984). As predicted, we measured no significant polarization in this case (see Fig. 9).
Spectral induced polarization measurements were performed using the apparatus developed at the Central Laboratory for Electronics, ZEL, at the Forschungszentrum Jülich, in Germany (Zimmermann et al. 2008). It is based on a four-electrode system operating in the frequency range 1 mHz-45 kHz. To our knowledge, Phase, [mrad] Test 1 Test 2 Figure 9. Test of the spectral induced polarization equipment and the electrodes (3M Red Dot) used to measure the potential response. This experiment is performed with a clean (clay-free) limestone sample saturated with a conductive brine (at 0.78 S m −1 at 25 • C). In this condition, we expect the phase to be close to zero in the investigated frequency range in agreement with the experimental data.
its accuracy is unmatched (0.1 mrad between 1 mHz and 1 kHz). We have tested the linearity of the system to check the influence of the input voltage on the measured currents (Fig. 10). This test shows that the spectral induced polarization mechanisms involved in these experiments are linear for applied voltage differences higher than 0.1 V. In all the experiment, we used U = 5 V for the applied voltage.
For each electrode configuration set-up, we need to convert the measured resistance into a resistivity by using a geometrical coefficient, which depends on the position of the electrode, their extension, and the geometry of the sample. Considering the sample and the measurement setup (Lg or Rd) discussed above, we determined the values of K by solving the Laplace equation with source current terms using the finite element code COMSOL Multiphysics 3.4 accounting for the geometry of the core sample and the geometry of the electrodes discussed above. The values of these coefficients are reported in Table 4.

Drying procedure
We dried the core samples in two phases: (1) We call the first phase the desaturation phase while (2) the second phase is termed the heating phase (see Cosenza et al. 2007 andGhorbani et al. 2009 for more details about this procedure). The studied samples were already partially unsaturated at the beginning of the experiments. We choose not to resaturate them to avoid possible perturbations like the swelling associated with changes in pore water composition. The saturations of the samples have been monitored through their weight with a precise balance (±0.1 g). The sample water saturations have been calculated from: where w = (m meas − m d )/m d is the gravimetric water content with m meas , the measured mass, and m d , the dry mass of the sample. Test 10 V Test 5 V Test 1 V Test 0.1 V Test 0.01 V ϕ Figure 10. Linearity test of the equipment used to perform the spectral induced polarization measurements in the frequency range 10 mHz-45 kHz. We used a partially saturated Callovo-Oxfordian core sample EST29300. (a) Amplitude of the resistivity versus the frequency. (b) Phase versus the frequency. We used the following input voltages: 10, 5, 1, 0.1 and 0.01 V. Note that the response is independent on the applied voltage (at least for voltages higher than 1 V). The responses conform to the linear Ohm's law. is the 'Total Combinable Magnetic Resonance' (TCMR) porosity, ρ w = 1000 kg m −3 , and ρ d = (1 − φ)ρ g are the water and the dry density, respectively (with ρ g = 2690 kg m −3 the grain density, see Gaucher et al. 2004 for details). We have decided to use the TCMR porosity (from borehole measurements) as a reliable estimate of the connected porosity because a study by ANDRA (2009) demonstrated that this was indeed the best estimate of the real porosity. The dry mass of the sample have been measured at the end of each experiment with the core samples being 24 hr at 105 • C in an oven. Table 4 presents the values of the porosity and the initial saturation, s i w , for each core sample.
We describe now quickly the desaturation phase. In this phase, the saturation decreases naturally by evaporation at T = 22 ± 2 • C and the relative humidity h r = 22 ± 4 per cent. Spectral induced polarization measurements have been done at different times after we removed the protective bags of the core samples, thus at different levels of the water saturations. The second drying phase is called the 'heating phase' because the samples have been put into an oven for 24 hr at three different temperatures (70, 90 and 105 • C). The higher temperature allowed the measurement of the dry mass of the sample, m d , for the saturation determination (see above). After each heating step, we have waited a couple of hours (typically 2-4 hr) before the spectral induced polarization measurements were performed. This is the time required for the sample to cool down at room temperature (T = 22 ± 2 • C). Fig. 11 shows the monitored saturation of the five samples during these two phases. Note the two phases for the desaturation: the desaturation phase and the heating phase.

I N T E R P R E TAT I O N
The spectral induced polarization signal can be decomposed into an amplitude |ρ * | and a phase ϕ. The behaviour of the amplitude can be compared to the DC resistivity except for the Stern contribution to surface conductivity (there is no contribution of the Stern layer at zero frequency according to the model described by Leroy et al. 2008 and). We first analysed below the amplitudes (at a given low frequency) and then later the amplitude and the phase versus the frequency.

Interpretation of the amplitudes
We first focus on the amplitude of the complex resistivity at 10 Hz.
We use the conductivity model developed by Revil et al. (2007) (eq. 15 above), which accounts for the saturation of the water phase. The conductivity of the pore water is calculated from the 220 D. Jougnot et al.
THERMOAR geochemical model (see Table 2) and the Donnan equilibrium model discussed by Revil & Linde (2006). In this approach, the contribution of the diffuse layer is accounted for inside the conductivity of the pore water phase. Using a 'non-DC' input current, we have to consider the Stern contribution to the surface electrical conductivity σ S . Therefore the conductivity model has five parameters: F, n, s c w , ϕ m and σ S ., which is itself determined from a grain size distribution. We determined the values of these five parameters for each core sample using the Simplex algorithm to solve a non-linear optimization procedure. This algorithm is able to optimize parameters value by minimizing a cost function G defined by the sum of the square residual between the N measured and calculated conductivities, σ k meas and σ k calc , respectively This problem being well posed, we do not need any regularization term. Miminization is performed for both the real and the imaginary parts of the complex conductivity. Fig. 12 shows that the Revil et al. (2007) model fit pretty well the experimental data with the optimized parameters given in Table 5. The optimized values of the different parameters agree with previous data obtained on core samples from the C2b1 and C2b2 units of the Callovo-Oxfordian formation. For example, F ∈ [38; 81] and ϕ m ∈ [ −49.5; −37.1]mV are consistent with the results obtained by Jougnot et al. (2009) and based on the diffusion of ionic tracers and in situ osmotic pressures in the Callovo-Oxfordian formation. The range of saturation exponents n ∈ [1.44; 2.04] is consistent with the values reported by  who developed a thermal conductivity model for unsaturated Callovo-Oxfordian clay-rocks. This model is an extension of the model developed by Revil (2000) to represent, in a unified framework, the electrical and thermal conductivities of granular media at saturation. The range of values of the saturation threshold s c w ∈ [0.10; 0.24] obtained from the conductivity measurements provide the first estimate of this parameter for the Callovo-Oxfordian clay-rock. The distribution of the surface conductivity associated with the Stern layer σ S ∈ [0.96; 3.25] × 10 −5 S m −1 is quite narrow. It is controlled by the density of counterions (sodium in the present case) in the Stern layer. As discussed further, this density is consistent with electrical triple layer calculations (see Leroy et al. 2007 and).
The two experimental setups (Lg and Rd) have allowed measurements perpendicular and parallel to the bedding planes (Table 5). From these considerations, we see that anisotropy has an influence upon n and s c w . In fact, the value of n normal to the bedding is higher than the value of n parallel to the bedding plane. The opposite occurs for s c w . However these observations should be considered with caution. Indeed in the Rd-setup, the current electrodes are very close to each other. This geometry implies that the current lines are not  homogeneously distributed through each core sample and then the measurements may not be representative of the complete volume of the core samples. The only way to check these observations would be to make electrical measurements with the Lg-setup on core sampled parallel to the bedding planes. This has been done by Cosenza et al. (2007) using the Tournemire clay-rock but not in this work using the Callovo-Oxfordian clay-rock.

Amplitude and phase versus frequency
We consider now the whole set of complex resistivity data and compare it to our mechanistic model for unsaturated clay-rocks. Fig. 13 presents the measured data for the spectral induced polarization measurements on the sample #EST29300 at seven different saturations. We can note the similarity in this trend with our model prediction as shown in Fig. 4. Indeed, as expected, the amplitude and the phase increase when the saturation of the water phase decreases. We also observed that after the heating phases, the experimental data are noisy and characterized by high uncertainties (especially after the heating step at 105 • C), which makes them difficult to interpret (Fig. 13). This could be explained by the increase of the contact resistance between the electrodes and the surface of the core samples.
In order to avoid misinterpretation of our results, we focus below only with the experiments performed with the Lg-setup excluding the data obtained in the heating phase (samples #EST27906, #EST29296 and #EST29300). In addition, in the heating phase, the formation of microcracks was observed and this creates a strong response upon the spectral induced polarization data Ghorbani et al. 2009) and the presence of these microcracks has not yet been accounted for in our model. In order to compare the experimental data and the proposed model for spectral induced polarization in unsaturated clay-rocks, we need to determine the unknown parameters of the model. In the previous section, we have determine four of the seven key parameters for the  five core samples: F, n, s c w and ϕ m (the values of these parameters are reported in Table 5). We used  measurements in saturated clay-rich porous media for the three other parameters ( 0 Na , τ 0 and c). These sets of values are considered to be prior values in the following. In order to find the best fit for each measurement, thus for the different saturations, we used a simple trial and error method from the prior parameters. Note that we optimize the whole complex signal: amplitude and phase at the same time.
The measurements on the sample #EST27906 (at s w = 0.82) have been fitted by the proposed model with the following parameters: F = 54.5, n = 1.79, s c w = 0.10, ϕ m = −44.1 mV, 0 Na = 0.891 nm −2 , τ 0 = 3.46 × 10 −2 s and c = 0.855. Fig. 14 shows that our model is able to fit pretty well the measured complex resistivities (total rms = 4.18 per cent with equal weight given to the in phase and quadrature conductivities). Dividing the frequency range in two parts, we see that the model perfectly fits the amplitude and the phase for frequencies higher than 1 Hz (Leroy et al. 2008;). At lower frequencies, the predicted amplitudes are still pretty close to the experimental data, but the modelled phase underpredicts the data. Fig. 15 shows the comparison between the measurements, during the desaturation procedure, and the prediction of the model. At low saturations, the measured Maxwell-Wagner polarization is higher than predicted by the model (from 1 to 15 kHz). The measured polarization is very well fitted between 1 Hz and 1 kHz.
Three explanations could explain the observed mismatch at low frequencies.
(1) Another mechanism (like membrane polarization, see discussion in ) may occur at low frequencies and has not been accounted for in our model. (2) A second possibility is that we have only considered sodium as the only counterion responsible for the polarization of the Stern layer. Maybe this assumption is too restrictive. Other cations are sorbed and may participate to the polarization of the Stern layer broadening the distribution of the relaxation times of the polarization of this layer in the frequency domain. (3) We have considered an unimodal distribution of the relaxation time distribution. However, the distribution of the relaxation times may be multimodal (e.g. Ghorbani et al. 2009 used a double Cole-Cole model to fit their data). The mean relaxation time, τ 0 , corresponds to a grain radius of 9.5 μm consistent both with the size of the tactoids and the size of the silica and carbonates grains.
The estimates of the model parameters determined by the best fits between the model and the data for core samples #EST27906, #EST29296 and #EST29300, are displayed in Fig. 16 as a function of the saturation. Note that the values of n, s c w and ϕ m optimized using the amplitude data at 10 Hz gave the best fit (Table 5). The optimized values of F and c are independent on the saturation (see Figs 16a and b, respectively). This is consistent with the fact that the formation factor characterizes the topology of the pore space (Revil & Cathles 1999) and the exponent c characterizes the broadness of the grain size distribution (Leroy et al. 2008). The mean values of the formation factor areF = 46.1 ± 1.8, 85.5 ± 5.1 and 44.1 ± 2.2 for the core samples #EST27906, #EST29296 and #EST29300, respectively. These values are consistent with those reported in Table 5. The mean values of the Cole-Cole exponents arec = 0.852 ± 0.005, 0.854 ± 0.009 and 0.827 ± 0.006, for the core samples #EST27906, #EST29296 and #EST29300, respectively.
We found that τ 0 and 0 Na vary slightly with the saturation. In fact, the clay mineral deformation (shrinkage) might break Stern layer connection inside the tactoids, and then increase the Stern layer diffusion path, thus the relaxation time. These results could confirm that the polarization of the Stern layer occurs in the tactoids and not around the coarse grains. This implies that the relaxation time distribution in such complex medium does not correspond directly to the particle size distribution, but more to a Stern diffusion path distribution inside the clay aggregates. In that case, τ would not be related to the grain radii but to the characteristic length of the tactoids. This explanation seems to be more realistic considering the complex geometry of clay assembly in consolidated clay-rocks (see Montes et al. 2004 for Scanning Electron Microscope pictures of the Callovo-Oxfordian at 10 and 50 μm scale). We observe also that the surface site density of sodium counterions in the Stern layer diminishes when the water saturation decreases. This could be due to the disconnection of the Stern layers inside the tactoid.

C O N C L U S I O N S
We proposed a mechanistic model of spectral induced polarization to unsaturated clay-rocks. Sensitivity tests were conducted to determine the influence of saturation and seven key parameters on the modulus and phase of the complex resistivity. The model was successfully tested against experimental data on partially saturated core samples from the Callovo-Oxfordian formation. The model agreed with the amplitude of the conductivity (at 10 Hz) as a function of saturation. Then, we showed that the model was able to reproduce the complex conductivity behaviour of the medium at high frequencies (>1 Hz) for water saturations comprised between 0.2 and 0.8. The polarization of the Stern layer that seems to take Na . Note that the optimized values of n, s c w and ϕ m are constant for all the saturations. The dependence of the relaxation time with the saturation is opposite to that observed by Binley et al. (2005). place at the scale of clay aggregates (tactoids) is related to a characteristic length that depends on the saturation of the water phase. In a future work, we will take into consideration the polarization of the Stern by several counterions with different valences and mobilities. In our opinion that would allow us to improve the fit at low frequencies (<1 Hz). We will also take into account the membrane polarization mechanism to produce a more complete theory.