Non-linear behaviour of electrical parameters in porous, water-saturated rocks: a model to predict pore size distribution

geo-resistivity meter. This paper also shows that the well-known Marshall and Madden’s equations based on Fick’s law cannot be solved without speciﬁc boundary conditions.


I N T RO D U C T I O N
Geophysical survey methods based on resistivity and induced polarization (IP) have been popular for a long time and applications in the subsurface, from ore and groundwater prospecting to geological research, pollution characterization, etc., abound. Resistivity measurements rely on the basic assumption that minerals and rocks obey Ohm's law. Likewise, IP was always considered, at least implicitly, as a linear phenomenon with chargeability (the parameter measured in the 'time-domain' version of the IP measurement technique) being independent of the excitation current. To the authors knowledge Bleil (1953) was the first to report on the non-linearity in IP measurements at variable current strength and charging time. He noticed this while studying laboratory scale models and, at that time, this phenomenon was known, discussed and accepted, although not properly understood. However since equivalent electrical circuit models have been widely used to describe the electrode polarization behaviour, which are implicitly 'linear', the two basic assumptions that have been accepted by nearly all geophysicists in the world are: (1) IP processes at time on and time off are the same. (2) There is a linear interdependence between applied electrical current and IP response. Practically all existing instruments record transient decay during time-off, that is in the absence of the main electrical field, when its measurement is easier.
However a new issue sometimes is a fully forgotten old issue.
results obtained in the laboratory on samples (Zadorozhnaya & Maré 2011). Moreover, since the model was strictly based on pore size distribution (PSD), it was able to estimate PSD on porous samples from geophysical measurements, as demonstrated by specific, direct measurement of pore size distribution of those samples, by means of a mercury porosity test. As is well known, pore size distribution of a porous rock is linked to its hydraulic conductivity (Schön 1996). Therefore, this finding is crucial for hydrogeological applications, among others, as it could offer the possibility to get estimates of the hydraulic conductivity by performing suitably relevant resistivity and induced polarization measurements. In this paper, the most meaningful laboratory data are reviewed and discussed, and then the theory of the predictive model is briefly recalled. A comparison between predicted pore size distribution of a sample and that measured by the mercury injection porosity test on the same sample is also presented. After this, a paragraph is dedicated to discuss why non-linear behaviour is of occurrence in field measurements that the geophysical community never had doubts that Ohm's law holds and never speculated about the nonlinearity phenomena in IP measurements. As the observation of non-linearity needs strong enough excitation currents and/or long charging times, it is relatively easy to observe this phenomenon, together with a violation of Ohm's law, on data acquired on samples in the lab, while in the field the ranges of supplied currents and charging times are, respectively too small and short, for practical reasons, to show any non-conventional phenomena.
Finally, bearing in mind that geophysical, non-invasive investigations are the most cost effective approach, a paragraph is dedicated to highlighting specific procedures for measurements at surface at a test site, where a sandstone aquifer is present and saturated with potable water. Purposely, a currently available commercial georesistivity-meter was used. The obtained results are discussed with a view of improving possible future measurement procedures during field practice.

P R I N C I P L E S O F T H E M E M B R A N E I P E F F E C T
The model for membrane induced polarization (diffusion coupling) was proposed by Marshall & Madden (1959). It was shown that the diffusion gradients and the electrical potential gradient had to be considered as the primary driving forces for ion motion. Anderson & Keller (1964) were the first to use and to apply the heat diffusion equation to explain the phenomenon of membrane polarization in rocks. It was found that the diffusion of ions through rock pores occurs as diffusion through a half limited capillary tube. This simplifies the equation to a homogeneous diffusion equation (at time-off). We can also refer to Schön (1996), Pipe et al. (1987) and Titov et al. (2002) who made contributions to the theory of membrane polarization as caused by the constrictivity of pores.
The foundation of membrane polarization is as follows: when electrical current flows through a channel containing pores with different radii (transfer numbers: see below eq. 2), an excess/loss of ions accumulates at the boundaries (Marshall & Madden 1959). If a pore space contains many parallel and negatively charged capillaries then the counter ions will be cations and the remainder ions will be anions. The cations moving to the cathode will pass the boundary between narrow capillary II and enter into the volume V A , at the anode side of the wider capillary I and then move further along (Fig. 1). Figure 1. (a) Model of a capillary: I and II are wide and narrow parts of the capillary, accordingly, l 1 and l 2 are their lengths, V C is the small volume of the capillary oriented towards the cathode, V A is the small volume of the capillary oriented towards the anode. (b) Distribution of salinity in the capillary after applying an electrical current (after Kobranova 1986).
Masses m IIk and m Ik (in moles) of cations and anions which enter into the volume V A and left this volume, are according to Faraday's law equal to (Fridrikhsberg 1995): where q is the amount of charge (in Coulomb, C) that passed through the boundary II-I, I is current, F is the Faraday constant (9.64955 × 10 4 C mol −1 ), z k and z a are valences of cations and anions, n I Ik , n I k are transfer numbers in narrow and wide capillaries. Subscripts k and a here and below indicate the cations and anions, respectively. Transfer numbers define the ratio of electrical current transferred by ions and can be described for cations and anions as: where M I k , M I Ik , M I a and M I Ia are effective (average) mobilities of cations and anions in the wide I and narrow II capillaries, respectively. Obviously cations are more mobile (transfer more electrical charges) when in a wide capillary because in narrow capillaries some of the anions are adsorbed by the double electric layers (DEL) and are immobile. In the narrower section, the transfer number of cations can be n I Ik → 1 and likewise n I Ia → 0. In the wide capillaries transfer numbers n I k and n I a identical and coincide with transfer numbers in free solution (in free solution n I k = 0.5 and n I a = 0.5). In this case the DEL does not play a considerable role for current flow. It was shown (Kobranova 1986) that if the surface areas of capillaries I and II are different and n I Ik > n I k then: This means that the salinity increases because cation concentration in the volume V 1A increases (see Fig. 1b). The concentration of anions is increased by the same amount (in case of a neutral solution). It was shown (Fridrikhsberg 1995) that at the boundary II and I, the salinity of cations and anions in volume V 1C decreases: Note that according to Kirchhoff ' law no charge accumulation occurs within an electrical circuit: It means that the amount of cations transferred through the contact pores with different surface areas is equal. Increase/loss in concentration of cations and anions levels off with distance from the contacts.

Mathematical description of membrane IP effect
The theoretical considerations about the distribution of ions in a solution filling the pores when a pulse of electrical current is applied were discussed by Zadorozhnaya et al. (2008Zadorozhnaya et al. ( , 2009Zadorozhnaya et al. ( , 2011. Let us briefly describe the model, equations and solution. The primary model consists of three connected pores with surface areas S 1 (central pore), S 2 and S 3 are left and right hand pores, respectively (Fig. 2).
The salinity distribution in a limited pore (central bar) is to be described using a diffusion equation. For the 1-D case (along x axis) the diffusion equation (Fick's Second law) is: where u is the salinity of the solution, a 2 = D is a diffusion coefficient, t is the time. Any problem in mathematical physics such as the above should be solved using specified boundary and initial conditions. The boundary and initial conditions for eq. (5) are: where l is the length of the central pore. Obviously in a stable condition (before electrical current will be switched on) ϕ on (x) = u 0 = const and equal to the initial salinity of cations and anions u 0 . It was shown by the above-mentioned authors that the boundary conditions can be written as follows: where Numerical subscripts indicate the number of pores (1, 2 and 3), u is the excess/loss of ion concentrations at the boundary between pores, σ is the conductance of ions, n is the transfer number of ions in the pore, I is the electrical current flowing in the pore (I k = I a ). The cause of the boundary conditions has been discussed in detail in (Zadorozhnaya & Hauger 2009). The solution of eq. (5) with boundary and initial conditions (6) can be found as a series (Koshlyakov et al. 1970): where Eqs (9) (boundary conditions) show a linear dependency of salinity on time at the boundaries between capillaries and depends on the square value of the electrical current I.
In accordance with the direction of current flow (Fig. 2) the ion concentration at the left side of the model increases while the concentration at the right hand side of the model decreases. Obviously a concentration decrease cannot continue infinitely, ion concentration at the right side of the model will reach zero causing a rupture of the electrical circuit. Then no electrical current flows through the capillaries any more, since the current pass is blocked. The conditions of the blocked capillary are: where t 0 is the time of blockage, which means the time when rupture of the electrical circuit occurred: numbers, conductivity of pore fluid and pore radii of the connected pores. The amplitude of the potential difference (voltage drop) also depends on the ion mobility and diffusion coefficient. The process of polarization continuous up to time t 0 , after which the rupture of the electrical circuit occurs and the potential difference between the ends of the pore becomes constant. During the polarization process all contacts between pores of different transfer numbers will be blocked and the electrical current will flow through the remaining pore channels. This brings us to the definition of the phenomenon of membrane polarization: Membrane IP is the successive blockage of interpore connections due to the excess /loss distribution of ions during current flow. Now we have to: (1) Demonstrate the non-linear IP effect in laboratory measurements.
(2) Prove the correctness of the proposed mechanism of the membrane induced polarization effect by a calibration with a mercury injection capillary pressure test (MICP).
(3) Theoretically and experimentally determine how the amplitude of the IP effect and resistivity depends on the amplitude of the applied current (current density).
(4) Show field results where geoelectrical parameters (i.e. resistivity and chargeability) depend on amplitude of current and duration of transmitting pulses.

Laboratory measurements
Membrane polarization occurs in all rock types if for connected pores, the surface areas and transfer numbers differ. L. Maré and M. Hauger studied thousands of samples of sedimentary rocks, ores, metamorphic and magmatic rocks and most of them demonstrated an IP effect. Some of the clastic sediment samples s (Dwyka age shale, shale with drop stones, mudstone and tillite), lava and samples containing non base metal mineralization (hematite and manganese ore) have been studied by Zadorozhnaya & Maré (2011). As an example the measurements on a sample of serpentinite metamorphic rock are presented here. Some samples of various sandstones will be analysed later in the paper.
The instrument RIP records voltage continuously due to current being injected into small rock samples as described below. The instrument was built by M.E. Hauger and has been used at the Council for Geoscience (CGS), Pretoria in their Physical Properties Laboratory where physical properties measurements have been performed for many years. The instrument contains a holder with two silicon medical electrodes (Fig. 3a). Several tests have been performed to prove that our electrodes are non-polarizable. A rock sample is located between the electrodes (Fig. 3b). At time 'on' both electrodes serve as transmitter electrodes, at time 'off '-as receiver electrodes.
The samples, small cubes 2.5 cm sided (or cylinders 2.5 cm in diameter), were dried in the oven for at least 24 hours and then immersed in tap water for several weeks. Then the measurements were performed using the two-electrode equipment for three to four different electrical currents from 0.0005 to 0.5 mA, depending on the rock's resistivity. The voltage drop was measured at the electrodes connected to the sample.
Since the geometrical shape of the sample is known, the resistivity can be calculated as follows: where V is the voltage drop, I is the supplied constant current, S is the surface area of the sample and L is the sample length. The chargeability is calculated as: where V 0.5 is the decay voltage measured at a time 0.5 s after cut-off and V p is the primary voltage. The method of measurements using the RIP instrument is described in detail in Zadorozhnaya & Maré (2011). In most of the cases a non-linear IP effect was present, however in some samples (soil containing animal manure, brown chrome, coal) only a linear dependence was observed. Fig. 4 shows the voltage and resistivity versus time and chargeability measured at different current strengths for a serpentinite sample, selected as an example of non-linearity. Minerals in this group are formed by the serpentinization of olivine rich rock, a hydration and metamorphic transformation of ultramafic rock. Serpentinite is formed from olivine via several reactions. Serpentinite does not contain quartz minerals, feldspar, mica, etc. The sample of serpentinite is a piece of a core collected in a borehole in one of the chrome deposit areas in Zimbabwe (sample C11068A). The magnetic susceptibility of this sample of serpentinite is 1.43 × 10 −3 SI, density is 3.415 × 10 3 kg m −3 and porosity 8 per cent. The sample shows a non-linear relationship between supplied signal and the response, both in terms of resistivity and chargeability. Note that at large time intervals the resistivity increases considerably with decreasing electrical current. However let's look at the details: at the earliest time there is a linear dependence of resistivity on current. This phenomenon will be discussed later. Chargeability always decreases with increasing electrical current (and current density).

Mathematical modelling of membrane polarization: calculating algorithm
If an electrical current I is flowing through the sample containing numerous pore configurations the total potential U can be presented as a sum of individual potentials due to the flowing of conductive current U curr (t) and excess/loss of concentration current U exc (t), unevenly distributed along the pores and related to the polarization effect in the rock: Then we have to calculate both parts of eq. (15). The model containing three connected capillaries, as depicted in Fig. 2, is too simple to be used for the interpretation of the laboratory measured data. A more complicated model is required. Let us consider a model containing parallel pore channels with m different pore sizes (Fig. 5). Then we will modify this model for a better approach to real pore structures.
The coded algorithm SAMPLE, calculating potential difference, resistivity of rocks and modelled pore size distribution in a sample has been written in Matlab. Eqs (12) and (13) contain numerous physical parameters, some of them are not constant, as they depend on temperature and pore radius and must be defined for each investigated sample. Physical parameters present in (9) that influence the amplitude of the membrane polarization effect are: resistivity of fluid, transfer numbers, zeta-potential, thickness of double electric layers (DEL), diffusion coefficient, dynamic viscosity of water, mobility of ions. Many of them vary depending on pore radii, temperature, water salinity, etc. Therefore these parameters have to be analysed one by one as a part of them is automatically calculated by the algorithm.
Input parameters to the computing code are also the size and porosity of the sample, the applied electrical current, the temperature in laboratory (usually kept at 25 • C). The samples are saturated by tap water with ion concentration 0.0047 mol l −1 (0.2632 g l −1 ). The number of pore sizes can be arbitrary; usually 70-100 and more than one size of pore can be used to determine the internal pore structure. The number of pores of the same size can also be selected arbitrarily.

Calculating voltage due to dc current
According to Ohm's law the potential that occurs in a model containing parallel or tortuous pores/channels (Kozeny-Carman model) is equal to: where N is the total number of pore channels in the surface area of the sample, n i is the number of pore channels of each size, and R i is the total resistance of each pore channel. Obviously the model containing parallel pores and channels is very far from a realistic model. This model describes the potential that occurs in tubes with a surface area that is equal to the average of the surface area of all channels filled by free water. A more realistic model of rocks can be assumed to be a fractal model. In this case the winding channels meet pores of different size. This fractal model was modified by Hallbauer-Zadorozhnaya (2013) following the Dakhnov's (1962) and Semenov's (1948) approaches. This model contains spheres of different radii coated by double electric layers (DEL). We specified the number of grains of different size. However there are double electrical layers that surround each grain. Obviously the coefficient of porosity ϕ will be determined by the number of fractions of spherical grains n. The resistivity of this model ρ m is equal to (Semenov 1948): where P P is a parameter of porosity (or structural factor) which is the ratio of the resistivity of the matrix to the resistivity of the material filling the pore spaces, ϕ 0 is an initial porosity, ρ f is the resistivity of the pore fluid. Our model is arranged by packing spheres of different diameters and irregular pores of different sizes into an arbitrary volume. Then the potential due to conductive current U curr (t) can be written as: Eqs (16) and (18) have been added to the code SAMPLE for calculating the potential difference due to the conductive electrical current. We have to note that during the time-on the U curr (t) will increase because of the increasing number of successfully blocked pores and decreasing numbers of non-blocked pores and then this has to be recalculated for every time where pore blockage occurs.
The same model (Fig. 5) has been used for calculating U exc (t) due to a salinity distribution in pores at each time sequence. This model is very effective because the distance, where excess/loss of concentration occurs, is much shorter than the length of the pore channels that connect the pores.

Resistivity of pore fluid
Due to the presence of the DEL effect the resistivity of the pore fluid is different from the resistivity of the confined electrolyte. The average conductivity of the pore fluid (including free solution and DEL) in polar coordinates (r, ϕ) with the reference located in the centre of the pore is equal to (Fridrikhsberg 1995): where where C 0 is the salinity of cations and anions in free solution (in mol l -1 ), C a and C k are the salinity of anions and cations (mol m −3 ), respectively, T is temperature ( o K) and R = 8.31467 is the gas constant (in J/mol·grad). The function 1 (r ) is the psi-potential, that is it is the potential in the inner plane of the DEL. The value of 1 (r ) is close to the zeta potential ζ . The 1 (r ) or ζ -potential must be specified a priori. How can the value of the ζ -potential be estimated? It was shown that the value of the psi-potential for ceramic diaphragms at the outer plane of the DEL 2 can reach ∼170-195 mV (Fridrikhsberg & Barkovsky 1984). Following Kormiltsev (1995) we regard the values 85-100 mV as optimal for clayey sediments and propose the following way for calculating the theoretical value of the ζ potential in the inner plane of the DEL.
Using (19) and (20) we calculated the conductivity of the electrolyte versus radii for each pore size and channel radii (Fig. 6). Fig. 6 shows that the resistivity of a pore filled with fluid and with a radius of less than 5 µm is significantly different from the resistivity of just a free electrolyte. For example if the thickness of the solid part of the DEL is equal to 3 ηm, the characteristic thickness of the DEL (for salinity 0.2632 g l −1 , ζ potential 0.1 V and temperature 18 • C) is 4.5 nm. Thus, the radius of 5 µm does not seem to be quite that narrow.
The fluid resistivity saturating the rock pores depends on the temperature. Usually it is assumed that: where ρ 18 • C and T 18 • C is the resistivity and the temperature at 18 • C, T• C is the calculated temperature. The SAMPLE code calculates the resistivity of the fluid for each size of pore and pore channel using the Figure 6. Average resistivity of pore fluid with temperature. Salinity is 0.0047 mol l −1 (0.263 g l -1 ). Curve index is pore radius (in mm). The resistivity of an unconfined electrolyte with the same salinity is also shown.
at Flinders University of South Australia on January 12, 2016 http://gji.oxfordjournals.org/ Downloaded from laboratory temperature. Underestimating the pore fluid resistivity can lead to a significant error in the mathematical simulation.

Calculating the overvoltage
For calculating the overvoltage U exc (t) each pore with unevenly distributed ions can be regarded as a single extended charged body.
That is why the model containing parallel pores and channels actually is very close to the real model. The potential U exp si po , created by a single pore l can be written as (Kalashnikov 1977): where ρ(t) is the volume of charge density (in C m −3 ), which is changing with time, r is the distance between the observation points and electrode (in m), ε 0 is the dielectrical permittivity of vacuum (in F m -1 ), ε is the relative dielectric permittivity (dimensionless), dv is a small unit of the pore volume, dv = S · l, S is the surface area of a pore, ρ (x) = Fzu (x). Basically membrane polarization is mathematically considered as the difference between the transfer numbers at each pore contact and connecting channel as described in eq. (9). The list of transfer numbers has been calculated in advance using a simplified Poisson-Boltzman equation described by Kormiltsev (1995). The diagram of calculated transfer numbers for cations and anions versus pore radius for different zeta potentials is presented by Zadorozhnaya & Maré (2011). To the SAMPLE code, the function 'transfer _numbers.txt' to define values of transfer numbers for each contact of pore and capillary has been added.
It is known that the diffusion coefficient in free water can be calculated using the Einstein-Stock equation: where v d is the dynamic viscosity of the fluid (in Ns m -2 ), r p is the radius of the particles (ions) (in m) and N A is Avogadro number (N A = 6.02322 × 10 26 k mol −1 ). The radius of 0.18 nm is assumed for a Na ion. However, this equation cannot be applied to fluid filled pores because the diffusion coefficient also depends on the pore radius r. Yaroshenko et al. (2004) presented an experimentally defined dependence D(r) as shown in Fig. 7: the diffusion coefficient decreases with pore radius. Usually channel radii are in the range from 14 nm to 5 µm. These authors also showed that the diffusion depends on the tortuosity (interconnectivity) and can be written as: where D 0 (r ) is the diffusion coefficient depending on the pore radii (from Fig. 7), T u is the tortuosity of pore channels (can be accepted as 1.5-2.5: see Yaroshenko et al. 2004), k v = v d tem /v d 20 , is the coefficient for the correction of water viscosity, v d tem and v d 20 = 1 are the water viscosity at laboratory temperature and at 20 • C, respectively. For each pore and channel D 0 (r ) is to be defined by using eq. (25). Water viscosity is also a parameter that depends on temperature (Fig. 8). In accordance with the temperature in the laboratory the correction value for water viscosity is introduced.
The mobility of Na ions is equal to 8.946e-8, and ion mobility of Cl = 4.695e-8 (in m 2 Vs -1 ) (as example, Bockris & Reddy 1998). In spite of using the parameter of mobility in our computations, this parameter does not influence the final result because there is a linear dependence of conductivity on mobility as given in eq. (20). Substituting (20) in (9) the mobility M will be removed. We note that M is a parameter of each type of ion but it is also constant for each kind of ion. However the rate of ion motion (velocity, V ion ) depends on the applied intensity of electrical current E, i.e. V ion = EM, (in m s -1 ) and differs between wide and narrow pores.
To obtain the total potential U exc we have to sum the potentials occurring in each single pore in the sample: where m describes the pore size and n m is the number of pores of size m.
In the SAMPLE code the boundary conditions and concentration distribution of ions along the pores are calculated using (9) and (10) for each size of pore. Moreover the radii of pore channels can be set differently in the algorithm: it means that one can select for instance 10 per cent of channels to have one size radius, 17 per cent of another radius, etc. Nonetheless if pore radii are too small and comparable with the radii of pore channels, then the influence of such pores cannot be considered: the process of polarization becomes very slow and the time of blockage tends to infinity. The structure of the code is quite complicated; it takes several seconds to compute one version of the model on a pc with a 2.2 GHz Intel Core i7 processor and 4 GB of RAM. Our experience shows that it is quite easy to compute a model which matches the laboratory measurements using only one current as shown below. However we have to calculate the voltage produced by different currents (at least three) to reduce equivalence.

M AT H E M AT I C A L M O D E L L I N G O F T H E I P E F F E C T I N S A M P L E S
Four samples, from Berea, Tennessee, Island Rust and Coconino, made available by Prof A.M. Binley, Lancaster University UK, have been measured in the petrophysical laboratory at CGS, where we performed the mathematical modelling of pore spaces using the above mentioned code called SAMPLE. These samples have also been studied by Baker (2001) who made the mercury injection capillary pressure (MICP) test and thus are considered by us as calibration models. The Berea, Tennessee, Island Rust and Coconino samples, are all sandstones. The results of the mathematical modelling of all four samples are shown in Figs 9(a)-(d). Three currents have been used for each sample. The pictures demonstrate a good agreement between mathematical modelling using the SAMPLE algorithm and laboratory measurements thereby proving the workability of our model.

Calculated pore size distribution and mercury injection capillary pressure (MICP) test
Characterizing the movement of water and contaminants in groundwater systems is challenging because of the limited ability to determine a high spatial density of measurements in a non-invasive manner. Since IP measures the electrical conduction and charge storage associated with the ion movement within the pore space of a rock, the technique provides a direct measure of the internal surface area of a porous medium and, we believe, a means of determining the pore size distribution (PSD) in a non-invasive manner. The challenge is to develop mechanistic models that link pore size characteristics to the IP response. We have to be able to calibrate and compare our results with the results of a method which directly determines the pore size distribution in samples. We have to examine the IP behaviour of a number of well characterised sandstones. We aim to build on this work by formulating a unified model for IP mechanisms, thus allowing indirect estimation of PSD in groundwater systems. The preliminary results of PSD modelling were presented by Hallbauer- Zadorozhnaya et al. (2010).
Capillary pressure can be measured using a variety of methods: centrifuge, nitrogen sorption, porous diagram, but the MICP method has the advantages of decreasing the measuring time, at Flinders University of South Australia on January 12, 2016 http://gji.oxfordjournals.org/ Downloaded from allowing a larger range of pore sizes to be handled, and the capacity to handle irregularly shaped samples (Sills et al. 1973). Following Baker (2001), we describe the mercury injection capillary pressure technique briefly. The standard sample is a three-inch long cylinder with a one-inch diameter. The samples were heated to a temperature of approximately 550 • C to drive out surface bound moisture and then dried at a constant temperature of approximately 150 • C. After placing the sample in a sealed chamber, the air was driven from the chamber. A tube containing mercury is connected to the chamber allowing one to measure the amount of mercury entering the chamber. Mercury was then pumped into the chamber, filling the entire chamber and a scale reading was made of the amount of mercury entering the chamber. To increase the pressure in the chamber, nitrogen is applied as a continuous pressure source. As mercury enters the evacuated pore space with increasing pressure, measurements are made of the lowering of the mercury meniscus as a function of pressure. A mercury hydrostatic head correction is also applied to the pressure reading. Mercury saturation is accomplished when either the sample chamber stops accepting mercury with increasing pressure or the injection pressure reached 59 400 psi [409.5486 MPa] (the maximum pressure able to be safely measured by the mercury injection apparatus).
The mercury injection curves of samples Berea, Tennessee Island Rust and Coconino, expressed as pore size with accumulated volume of total porosity are presented in the Figs 10-13.
The porosity of the Berea sand-stone sample is 19 per cent with a mean grain size of 0.3 mm. Fig. 10 demonstrate the results of modelling of pore space for the Berea sample: cumulative pores and histogram of PSD. The agreement of cumulative pore curves obtained by MCIP and mathematical modelling is very good (Fig. 10a). The histogram of pore size distribution is shown in Fig. 10b. The pore size varies from 14.2 to 2.6 µm. However the pore sizes available for modelling are bounded from both sides. The widest pores are blocked at the earliest times. Then there is uncertainty: all wide pores block almost simultaneously. On the other hand the narrow pores are blocked at times reaching tens of minutes or hours. The range of pore sizes giving individual contributions to the mathematical modelling varies from tens to single micrometers, usually a span of 1.5-2 decades. In Figs 10(b)-13(b) the vertical axis corresponds to the ratio of simulated pores to the total volume of pore (porosity), that is V porei /V total pore . For a comparison of the shape of PSD for different samples we have plotted them on the same coordinate system. For clarity the curves are represented on a logarithmic scale in the diagram. Also for clarity we presented all curves as separated series. Returning back to the Berea sample we can see that the pores with a radii of 6-9 µm. dominate in this sample.
The porosity of the Tennessee sample is 6 per cent, the main grain size is 0.2 mm without interconnected pore space. Results of the mathematical modelling in comparison with the MCIP test are presented in the Figs 11(a), (b) and in the Fig. 14. The pore size distribution for large pores has an arched shape; it means the number of wide pores decreases sharply. The average of pore sizes shifted towards the narrower pores. The agreement of curves relevant to cumulative pores obtained by mercury injection capillary pressure test and mathematical modelling is again very good. Fig. 14 shows that wide pores occupy only a small volume of pore space. Most of the pores are less than 5 µm and cannot participate in the mathematical simulation of pore space.
The porosity of the Island Rust sample is 14 per cent. The main grain size is 0.2 mm with pores that are not well connected. However the MCIP and mathematical modelling curves show the best agreement of all four samples (Fig. 12a). The histogram of pore size distribution is shown in Fig. 12(b). The sample contains a large range of pore sizes and reaches 35 µm (Fig. 14). However small pores constitute a significant portion of the pore spaces of the sample.
The porosity of the Coconino sample is 11 per cent, the main grain size is 0.5 mm (Baker 2001). The Coconino sample is the only one which defied our modelling. The cumulative curve obtained by modelling has shifted towards the wider pores. We can only assume that the presence of iron oxide could affect the process of overvoltage relaxation. According to the modelling, the pores in the sample are large (Figs 13b and 14). The mean pore size is about 10 µm.

D E P E N D E N C Y O F R E S I S T I V I T Y V S C U R R E N T D E N S I T Y: D I R E C T O R R E V E R S E ? O H M ' S L AW VA L I D I T Y A N D L I N E A R I T Y O F T H E I P P H E N O M E N O N I N F I E L D M E A S U R E M E N T S : D I S C U S S I O N
We should expect good similarities between laboratory measurements and field experiments, that is we should expect to obtain data proving that chargeability depends on applied electrical current, wherever aquifers with internal pore structures and textures as described above and corresponding to the studied samples are present in the subsurface. However the field experiments show that resistivity and chargeability practically do not depend on the current. Rather this behaviour follows the universally accepted basis of the whole domain of geo-electrical methods, and upon it the repeatability principle of the measurements is established. It is worthy to remember that repeatability of measurements means that we must obtain the same results under the same measuring conditions: using the same electrical current is one of the measuring conditions. So how can we explain the dichotomy between the laboratory and theoretical calculations on the one hand that match fairly well and the field experiments on the other hand? The reason can be the following. Many years ago, when the IP method was rapidly growing, the geophysical community soon became aware that a non-linear relation between chargeability and excitation current could occur (Bleil 1953). This phenomenon was considered a 'saturation' effect, without any more speculation. Still three decades later, one of the current authors contributed to a paper (Iliceto et al. 1982), where the authors declared that their laboratory measurements on many samples of soils (clay, silt and sand) were all done avoiding non-linearity by means of a careful limitation of the injected current in each sample. Now we have to explain this contradiction between the laboratory measurements and our theoretical model with Ohm's law. As we mentioned previously the RIP instrument allows measurements over a wide range of electrical currents from 0.05 µA to 1 mA. We selected an appropriate sample (for example sandstone) for measuring at different currents (CP9012AA). Fig. 15 shows the resistivity dependence on current density (current divided by the surface area of the sample, A m -2 ). The numbers on the curves indicate the time of relaxation (from 0.2 to 3 s). The figure clearly shows that the resistivity depends on the current density, violating Ohm's law. At low current density resistivity linearly increases with increasing current density; moreover increasing the time of relaxation, this dependence becomes stronger. All curves reach a maximum, and then we observe a decrease in resistivity with increasing current density. Obviously if the current density is small then the accumulation of ions occurs very slowly and the IP effect is negligible compared to the dc voltage V DC [V I P V dc (I )]. In this case Ohm's law holds. If the current density increases then the contribution of IP also increases (V I P (t) < V dc (I )) and resistivity increases considerably. However, according to the definition of the membrane IP effect (Zadorozhnaya & Hauger 2009), V I P (t) tends to be constant with the increasing time of the supplying current ( lim t→∞ V I P (t) → const) and the time of blockage of the pores directly depends on the applied current (in square, I 2 ). It means that, at sufficiently high current densities, the ratio V I P (const)/V DC (I ) will decrease and the resistivity of the sample will also decrease as shown in the Figs 15 and 16(a) and (b). These results completely correspond to the proposed theoretical model of membrane polarization.
It is not practical to measure with such a wide range of currents all rock samples using the RIP instrument. However we can execute a mathematical simulation using the above described calculating  (b) show the theoretical dependence of resistivity on current density calculated for samples Tennessee (a) and Coconino (b). As we see, the curves have the same shape: a direct dependence of the resistivity for small current density and an inverse relationship at higher current density.
In the laboratory experiments and theoretical calculations illustrated above, the electrical currents range from 0.5 µA (the smallest current which ensured an acceptable signal-to-noise ratio) to 1 mA. We should recalculate this range in terms of current density to make these measurements comparable with field practice. For samples consisting of small cubes of 2.5 cm diameter size the current density covered the range 8 × 10 −3 to 1.6 A m −2 .
For the sake of comparison with field practice, we now calculate the current density distribution in the subsurface generated by a current dipole, made-up of two electrodes (A, B) placed onto a homogeneous half-space model of conductivity σ . The current density J is easily calculated at each point P from the potential U(P) of the dipole.
where I is the injected current). Based on Ohm's law: Current density J is independent of conductivity σ , which is obtained when combining eqs (26) and (27). In Fig. 17 the results of such a simulation on a vertical section below the current dipole is shown, for a dipole length of 20 m. With the exception of the pixels nearest to the electrodes, the current density is generally much less than 0.01 A m −2 (green to blue colour), which is of the order of magnitude of the lowest values used in the laboratory survey. The currently available multi-electrode geo-resistivity meters inject currents not greater than few A into the subsurface because of safety reasons and portability. Thus in field surveys the geophysical community observes that both Ohm's law holds without exceptions and non-linear IP is not experienced.
Chargeability values due to membrane polarization are usually small and do not exceed several mV V -1 . Fig. 18 shows the calculated signals (normalized) obtained for different current densities. This result shows that the overvoltage, calculated using (26) at low current densities, that is under usual field conditions, will always be very small and more or less proportional to the applied current (i.e. the 'linear' behaviour). The red ellipse indicates the segment of low current density and of short on-time ranges where the IP data are usually recorded in the field. Moreover these results confirm that we cannot expect valuable information on pore size, etc. to be extracted from chargeability whenever the current density is low. High current density and/or much longer charging time are required.
Thus using the laboratory measured data and mathematical modelling we eliminated the contradiction between the proposed model of membrane polarization and field experiments.

F I E L D E X P E R I M E N T
Applying to field practice the results described above should allow for a remarkable progress in hydrogeophysics, since characterization of possible non-linearities in electrical resistivity tomography (ERT) and IP measurements in the field should allow an estimation also of the PSD, besides the effective porosity, and thus, in an intermediate way, of the hydraulic conductivity of the aquifers present in the subsurface.
Nevertheless, before anything else we should briefly examine the existing differences between data acquisition in the lab and in  the field. Lab data were obtained, as presented above, by injecting selected and constant current intensities or, equivalently, current densities, so that each sample was characterized and its own PSD could then be estimated from the model. In the field, volumes of aquifer bodies are crossed by current densities which depend on many factors: injected total current, depth of the bodies, ratio between their resistivities and those of surrounding bodies and, last but not least, reciprocal positions of injecting current electrodes at surface and discernible bodies. Therefore the current density distribution cannot be directly controlled and the hypothesis to repeat acquisition using different injected currents, in compliance with the lab procedure, should therefore be carefully speculated on. Ideally, for each volume element (pixel or voxel) of the estimated resistivity distribution, we should know the current density that crossed it, for each current injection. Of course this datum can be readily established using eq. (27) as the resistivity distribution, that is the estimated model, is known and the electric potential drop on each pixel/voxel can be calculated in the whole volume by the modelling/inversion code for each position of the current electrodes. In this way we should recover the required information, but unfortunately it is not made available by current commercial inversion codes, nor would it be actually meaningful, as it would be derived by a distribution of resistivity values which was defined implicitly as 'current-independent' by the embedded modelling algorithms.
However, before analysing further this relevant topic, it must be shown that non-linear behaviour is an actual occurrence, since in all the geophysical literature linearity is presumed and no similar occurrences were reported up to now. Therefore the field experiment we present in the following is devoted to show that such an occurrence is realistic. It was planned in full analogy with our lab experience. We wanted to obtain localized information from a behaviour in the above described sense, by varying the injected current between two subsequent acquisitions, we obtain a significant variation of current density on the whole investigated subsurface and thus on each pixel, so that if a dependence exists, it should be observed by comparing the two resistivity/IP models. Moreover, bearing in mind the data shown in Fig. 18, we selected a time cycle of 1 s on, 1 s off, which is a common choice for time domain IP data acquisition.
A readily accessible site was selected, among the many sandstone formations present in the Northern Apennines (Northern Italy), where potable water is being exploited by an aqueduct plant. The obtained results are illustrated in Figs 20 and 21 and discussed in the following paragraphs. From a geological viewpoint the investigated area is part of the Emilia-Romagna Apennine (ERA) which constitutes the northern Apenninic rim of the Padan-Adriatic margin. It forms a complex sequence of fold and thrust belts resulting from the collision of the Adria Plate with the European Plate that started in the Late Cretaceous. Exposed rocks belong to Ligurian Units and in particular to the upper part of the External Ligurids. These include numerous Late Cretaceous and Palaeogene tectonic Units. The local outcropping flysch of Monte Cassino (mcs), Upp. Campanian-Maastrichtian (Papani & Zanzucchi 1969), is partially covered by a thin layer of detrital sediments belonging to ancient landslides (a3) and slope sediments (a2g, Fig. 19). The mcs unit is composed of a thick, monotone marly-calcareous, turbiditic sequence (Helminthoid Flysch Auct.) and of a matrix supported or clast supported detrital deposits (Basal Complexes Auct.). In the studied area, the prevailing outcropping rocks are composed of medium to fine grained sandstone layers whose base is locally overlaying fine grained sediments (clay).
Data were collected using a geo-resistivity meter model SAS4000, manufactured by ABEM Instruments AB (Sweden). This equipment is particularly well suited to this test since it is a 'current at Flinders University of South Australia on January 12, 2016 http://gji.oxfordjournals.org/ Downloaded from generator', that is the equipment injects the current intensity as selected by the operator throughout all measurements. Two profiles for carrying out electrical resistivity and IP measurements in a tomographic way were laid out and driving into the ground 25 electrodes spaced 3 m apart and 26 electrodes spaced 2 m apart, respectively. Taking into account that PVC pipes are used for collecting water and wells were build by bricks, no passive noise could affect the data and Profile 1 was purposely laid out very near to some of the water wells (solid squares in Fig. 19), so that we were reasonably sure that the aquifer was present in the area covered by our two profiles. The Wenner-Schlumberger array was selected, since this is the best combination of signal-to-noise ratio and of sensitivity to lateral discontinuities. In the time interval of current off, five adjacent integration windows were selected to measure apparent chargeability M a , defined as: After waiting 10 ms after current shut-down, lasting 20, 40, 80, 160 and 320 ms, respectively M a is then measured in mV.s V -1 , that is in ms. Here we will present the data pertaining to the earliest window. Two measuring cycles were then carried out on each profile, using 50 and 500 mA, respectively, without changing any electrode position. Acquisition of the whole data sets together with topography required a day of field work for three people. Two tomographic sections of electrical resistivity and IP were then obtained after data inversion by using commercial software RES2DINV (Loke 2011), based on the Gauss-Newton method (Loke & Barker 1996).
Here we will show the results for Profile 1, directed SE-NW, because the results of the shorter Profile 2 were less significant, as it intersects Profile 1 where the aquifer is absent and on its western part it lies above the loose sediments of the local valley.
In Fig. 20 the resistivity section is shown, where the low-current results are reported in the upper window and the high-current results below, respectively. The inversion process converged in three iterations, with little difference in per cent rms errors between the low current case (3.5 per cent) and high current case (3.0 per cent). A few isolated data points, whose per cent difference with the predicted values were more than 10 per cent were excluded using the trim option of the software.
As we must transform this geophysical model into a geologicalhydrogeological one, we should remember both the geological information from literature, synthetized above. At this site two types of sediments outcrop: clay and sandstone. The latter appeared as coherent rock in wells and as stones of various sizes, dispersed at the soil covered surface or mixed with a clay matrix. The clay, saturated with fresh water, is generally a good conductor, with resistivities commonly spanning the range 5-20 Ohmm (see any textbook on Applied Geophysics), while the sandstone is much more resistive, even if saturated with (fresh) water. As a consequence the resistivity section can be transformed into a local geological/hydrogeological model where (i) areas coloured light to dark blue correspond to clayey sediments, (ii) the areas coloured in red and marked with (A) correspond to the aquifer (saturated sandstone), and iii) the mixture of clay and sandstone areas corresponds to green-yellow colours. It is likely that, in the most elevated portion of the profile, areas with a deep red colour, marked with (B) corresponds to arenaceous, drained slope debris, which should explain why these areas do not show any variation in resistivity with changes in the injected current mode. It is evident at a glance that areas marked with an A in the lower section of Fig. 20 are characterized by greater values in resistivity when compared to the upper section. First of all, it may be worth to estimate at least the order of magnitude of current density inside the body marked with (A) around an abscissa value (distance along profile) of 30 m, at a depth of about 10 m below surface (b.s.l.) by using the same procedure which led to results shown in Fig. 17. Taking into account that the depth of investigation of the Schlumberger array is about 0.2 times the distance between current electrodes A and B (in an homogeneous half-space) we repeat the calculations for A-B distance of 50 m for both injected currents of 50 and 500 mA. Then we obtain an average value of current density at 10 m depth of 0.00005 and 0.0005 A m −2 , respectively. These amounts correspond to lab measurements carried out with currents of 0.00003 and 0.0003 mA, respectively on samples of cross-sectional area of 6.25 cm 2 , but 500 mA is the smallest current density used in the lab measurement approach.
To highlight volumes where the resistivity varied, we constructed a section depicting a per cent variation calculated as (resistivity at 500 mA -resistivity at 50 mA)/resistivity at 50 mA), which is shown in Fig. 21. The strongest increase corresponds to the areas marked with an A in Fig. 20 and exceeds 20 per cent. A decrease of resistivity with increasing current is also noticeable, which corresponds to the resistivity volumes associated with clay. This effect should be investigated further apart, as it is not predicted by the above-discussed model, which concerns only porous textures.
Similarly in Fig. 22 we show the IP tomographic sections. As can be seen in Fig. 22 the decrease of chargeability with increasing current is dramatic and generalized to the whole section. Therefore the effect of reducing chargeability with increasing current occurs not only within the aquifer bodies but also, at least, within the mixed clay-sandstone volumes. Before any further investigation of this ubiquitous difference can be performed, perhaps a contribution from ambient noise should be taken into account, considering the more unfavourable signal-to-noise ratio of the low current data, which could be, at least partly, responsible of the general higher values, while the red-brown maximum in the low current section is preserved, although strongly weakened, in the high current section.

D I S C U S S I O N A N D C O N C L U S I O N
In this paper, two separate but related goals are tackled. The first one is to demonstrate that, in some specific saturated rock textures, non-linear behaviour of IP and violation of Ohm's law not only are a real phenomenon, but they can also be satisfactorily predicted at Flinders University of South Australia on January 12, 2016 http://gji.oxfordjournals.org/ Downloaded from by a suitable physical-mathematical model. This model is based on Fick's second law as the Authors showed that the well-known Marshall and Madden's eqs (37) and (38) (1959) are only correct to a certain point (see the Appendix for more details). Many tests on laboratory samples show the correctness of this approach.
The second is that as the model links the specific dependence of resistivity and chargeability of a sample of the injected current to its PSD, it is able to predict PSD from measurements that are current and time dependent and are in good agreement with MICP test results.
These facts open the perspective for hydrogeophysical applications in the field. To show the possible portability of the above results from the laboratory to field practice using ERT and IP methods, a specific field test was conducted. A suitable site was selected, where the subsurface hosted at shallow depths a sandstone aquifer, as shown by knowledge of the local geology and by the presence of several wells for potable water extraction. The measuring strategy was planned bearing in mind that the most important information should be available in a relatively short time. Evidence of the possible dependence of resistivity and chargeability on injected current was obtained, by acquiring two data sets per profile at two injected current strengths and resulted in a one order of magnitude difference chargeabilities, keeping the cycle time constant. A commercially available geo-resistivity meter was purposely used, to verify that the expected behaviour could be acquired using current technology. The obtained results can be summarized as follows, which are in excellent agreement with the proposed model: 1. A fair increase in estimated resistivity with increasing current within those volumes that can reasonably be associated to be saturated sandstone bodies.
2. A strong decrease in estimated chargeability with increasing current, which extends also over other volumes (besides sandstone) of the subsurface; this last finding should be investigated further.
In the interest of establishing the above phenomena for a larger body of case histories and experiences the authors will make data available to anybody interested for further verifying these findings. They are acutely aware that field evidence should be strengthened both by investigating other potentially suitable sites (there are many available in the Northern Apennines, similar to the one investigated here) and by investigating the effects of a longer current on and off time cycles, although this practice could be more time consuming.
Finally, we want to highlight that the present paper is not intended as a critique on the well consolidated geoelectrical and IP methods, which would be a misunderstood consequence of the proof that non-linear phenomena can take place. Rather it is the exploration of new, and potentially powerful features of these methods in hydrogeophysics such as, acquiring data (i.e. apparent resistivity and apparent chargeability) with commonly available geo-resistivity meters in a suitable way to highlight possible non-linearity, where quantitative information about PSD can be collected, which extends the characterisation of hydraulic properties of an aquifer well beyond the estimation of its effective porosity.