Seismoelectric response of heavy oil reservoirs: theory and numerical modelling

SUMMARY
We consider a linear poroelastic material filled with a linear viscoelastic solvent like wet heavy oils (oil as opposed to water or brines is wetting the surfaces of the pores). We extend the electrokinetic theory in the frequency domain accounting for the relaxation effects associated with resonance of the viscoelastic fluid. The fluid is described by a generalized Maxwell rheology with a distribution of relaxation times given by a Cole–Cole distribution. We use the assumption that the charges of the diffuse layer are uniformly distributed in the pore space (Donnan model). The macroscopic constitutive equations of transport for the seepage velocity and the current density have the form of coupled Darcy and Ohm equations with frequency-dependent material properties. These equations are combined with an extended Frenkel–Biot model describing the deformation of the poroelastic material filled with the viscoelastic fluid. In the mechanical constitutive equations, the effective shear modulus is frequency dependent. An amplification of the seismoelectric conversion is expected in the frequency band where resonance of the generalized Maxwell fluid occurs. The seismic and seismoelectric equations are modelled using a finite element code with PML boundary conditions. We found that the DC-value of the streaming potential coupling coefficient is also very high. These results have applications regarding the development of new non-intrusive methods to characterize shallow heavy oil reservoirs in tar sands and DNAPL contaminant plumes in shallow aquifers.


I N T RO D U C T I O N
In the last 70 yr, a number of laboratory, field, and theoretical investigations have established seismoelectric (seismic-to-electric) and electroseismic (electric-to-seismic) conversions associated with the charged nature of natural and synthetic porous materials filled with a solvent (usually a Newtonian electrolyte). Preliminary investigations have attempted to image the subsurface and to locate potential oil and gas reservoirs and ore bodies using these properties (Martner & Spark 1959;Migunov & Kokorev 1977;Kepic et al. 1995;Russell et al. 1997). These methods were also used more recently in hydrogeophysics to detect wet ice in glaciers (Kulessa et al. 2006), to localize fractures in a fractured aquifer (Fourie 2003), to image interfaces within the vadose zone (Dupuis et al. 2007), and to monitor electromagnetic disturbances associated with earthquakes (see Nagao et al. 2002;Huang 2002;Huang & Liu 2006). Measurements can be performed not only at the ground surface of the Earth 782 A. Revil and A. Jardani numerous examples in the literature where oil is the wetting phase for the solid grains (see Anderson 1986;Buckley & Liu 1998, for some field examples). We are not aware of any work in this direction so far. At low frequencies, it has been observed that electrokinetic phenomena in oil-wet porous media are similar to those existing in water-saturated porous materials (Yasufuku et al. 1977;Alkafeef et al. 2001;Alkafeef & Smith 2005). Various authors have shown that the components of oil that are responsible for wettability are also polar components (Buckley & Liu 1998;Alkafeef et al. 2006). An oil with such polar molecules is also a good solvent (Delgado et al. 2007). We will show that the theory developed in this paper is also consistent with recent developments in the theory of streaming current/streaming potential associated with the flow of the pore water in a deformable porous material (see Revil & Linde 2006;Crespy et al. 2008).

Properties of the fluid phase
We start this section by a description of the constitutive model of the stress-strain relationships for a viscoelastic fluid such as a wetoil. A linear Maxwell fluid model consists of a linear dashpot in series with a linear spring. In this situation, the fluid behaves like a Newtonian viscous fluid at low frequency (or for long timescales) and like a solid at high frequencies (or short timescales) (Behura et al. 2007). The transition frequency is discussed below. We note T i j and e i j the components of the stress and deformation tensors T f and e f of this fluid (subscript f ), respectively. Both T f and e f are symmetric second-order tensors. The stress components of the elastic and viscous contributions of the Maxwell fluid are given by where the superscripts 'e' and 'v' stands for 'elastic' and 'viscous' contributions to the stress tensor and to the deformation tensor, the dot above the symbol denotes the first time derivative of the property, λ f and G f are the Lamé and shear moduli of the fluid, η f is the dynamic (shear) viscosity and p f is the local fluid pressure. The bulk modulus of the fluid is defined by the relationship − p f I = K f (∇ · u f )I where I is the identity tensor and u f is the displacement of the fluid phase. It is related to the Lamé constants by K f = λ f +(2/3)G f . For a Maxwell body, the components of the stress tensor of the fluid, T f , are given by T i j = (T i j ) e = (T i j ) v anḋ e i j = (ė i j ) e + (ė i j ) v . Using eqs (1) and (2) Using harmonic variations of the stress (the stress oscillates in time as e −iωt ), we can take the previous equation in the frequency domain using a Fourier transform (we keep however the same notations). The characteristic frequency is defined as ω m = G f /η f and the associated relaxation time is τ m = η f /G f . Typically, for the Mexican crude oil investigated by Dante et al. (2007), the critical frequency depends on the temperature and is in the range 10-100 Hz in the temperature range 20-40 • C. For heavy oils in tar sands for example, the critical frequency of resonance can occurs at much lower frequencies. Usually, heavy oils are heated to decrease their dynamic viscosity and the resonance frequency can occurs in the seismic frequency band (see Behura et al. 2007).
In the frequency domain, eq. (3) is written as, At low frequencies ω ω m so that p * f = p f , η * f = η f , the third term on the right-hand side of eq. (4) is negligible and the stress tensor is given by (e.g. De Groot & Mazur 1984) where the (ḋ f ) v is the viscous contribution to the fluid strain deviator (e.g. De Groot & Mazur 1984), where d f is the fluid strain deviator, v f denotes the velocity of the fluid and where the superscript T means transpose. The deformation tensor is related to the deviator by (e.g. De Groot & Mazur 1984), At low frequencies the fluid behaves like a Newtonian fluid. At high frequencies ω ω m , the stress tensor is given by and therefore at high frequencies the pore fluid behaves as a solid. A perfect Maxwell fluid has a relaxation time distribution described by a delta function. Such a model can be generalized using a distribution of relaxation times, for example a Cole-Cole distribution [see for example,  for deformation of rocks by pressure solution, Cosenza et al. (2007) for spectral induced polarization]. Behura et al. (2007) used a Cole-Cole model to fit the hydromechanical behaviour of the Uvalde heavy-oil. A comparison between the Maxwell and Cole-Cole behaviour of the resonances of a generalized Maxwell fluid is shown in Fig. 1 where we compare experimental data from Castrejon-Pita et al. (2003) for a Maxwell fluid (a mixture of cetylpyridinium chloride and sodium salicylate) with a model in which the Cole-Cole distribution of relaxation times is used. From this figure, it can be conclude that the Cole-Cole model is able to represent the hydrodynamic behaviour of the fluid quite well. The probability distribution (normalized) of the relaxation time for a Cole-Cole distribution is given by (Cole & Cole 1941) In a semi-log plot, this distribution is symmetric with respect to τ = τ m , which corresponds to the peak of the distribution. For 0.5 ≤ c ≤ 1, the Cole-Cole distribution of relaxation times looks like a log normal distribution The tail of the Cole-Cole distribution is however increasingly broad as c decreases. This type of distribution has been used recently by Leroy et al. (2008) and Leroy & Revil (2009) to model the grain size distribution of packs of glass beads and sedimentary rocks to compute their complex resistivity.
The behaviour of the Maxwell fluid can be considered from the standpoint of the behaviour of an equivalent Newtonian viscous Figure 1. Hydrodynamic response of a viscoelastic fluid in a linear tube of radius R = 2.5 cm submitted to harmonic fluid pressure variations. The fluid is a mixture of cetylpyridinium chloride and sodium salicylate (CPyCl/NaSal 60:100) (data from Castrejon-Pita et al. 2003). The experimental data are scaled with respect to the viscous relaxation time τ = R 2 ρ f /η f (η f = 60 Pa s, ρ f = 1050 kg m −3 ).ξ denotes the normalized conductance (relative to the DC-value) and J 0 is the Bessel function. (a) Comparison between the experimental data and the response predicted by a linear Maxwell fluid (τ m = 1.9 s). (b) Comparison between the experimental data and the response predicted by a generalized Maxwell fluid using a Cole-Cole distribution of relaxation times (τ m = 2.9 s, Cole-Cole exponent c = 0.9, parameters optimized with a Newton-Raphson minimization method using the L 2 -norm). fluid with frequency dependent bulk modulus and dynamic shear viscosity (Behura et al. 2007) Using Boussinesq approximation, the flow of a fluid inside a deformable porous material is governed locally by the Navier-Stokes respectively, where ρ f denote the local mass density of the fluid, F f is the body force, and represents the viscous (deviatoric) stress tensor of the fluid. For a purely viscous Newtonian fluid, we have ∇ · π = η f ∇ 2 v f . Assuming a linearized Maxwell model and using eq. (4), the viscous stress tensor π is given bẏ In the frequency domain, eq. (21) yields Inserting eq. (23) inside the Navier-Stokes equation, eq. (17), yields Consequently, we come to the important conclusion that the Navier-Stokes equation of a viscous Newtonian fluid can be used by replacing the classical viscosity η f by an effective (time or frequency dependent) viscosity η * f . This is the main result of this section with the new result of using a Cole-Cole distribution of relaxation times to describe a generalized Maxwell fluid (see Fig. 1b).

Properties of the solid phase
We consider that the solid phase is formed by a monomineralic isotropic solid material assumed to be perfectly elastic. The local elastic equation of motion for the solid phase is where u s is the displacement of the solid phase, ρ s is the bulk density of the solid phase and F s is the body force applied to the grains. The microscopic solid stress tensor is given by where d s is the solid strain deviator, K s is the bulk modulus of the solid phase assumed to be isotropic and G s is the shear modulus.

P RO P E RT I E S O F T H E P O RO U S M AT E R I A L
We derived below the coupled constitutive equations between the Darcy velocity (Section 3.1) and the total current density (Section 3.2) for a linear poroelastic material saturated by a polar (wet) oil ( Fig. 2). We look for an average force balance equation on the fluid in relative motion with respect to the solid phase. We introduce the relative flow velocity v = v f −v s wherev s =u s is the phase averaged velocity of the solid phase. In the frequency domain, using eq. (24) yields where c is the Cole-Cole distribution corresponding to the distribution of the relaxation times for the generalized Maxwell fluid. The body force applied to the pore fluid corresponds to the electrostatic force associated with the existence of the local net charge density Q (expressed in C m −3 ) plus the body force due to the gravity field where ρ f is the mass density of the pore water, g the acceleration of the gravity and E is the electric field. The average forceF f can be also associated with a source mechanism (see below Section 6.4 and using the decomposition provided by eq. 56). We average now eq. (28) over the fluid phase. This yields

785
The flow is also subjected to the no-slip boundary conditionv = 0 at the surface of the pores.
The local charge density Q depends on the position inside the pore space. This local charge density is related to the concentrations of the ionic species that are affected by the Coulombic field created by the effective surface charge density of the mineral surface. These surface concentrations are described by Boltzmann distributions. The electrical double layer is always thin with respect to the size of the pores and the size of the grains. In this paper, we use the approximation developed by Revil & Linde (2006) that Q can be replaced by its phase average over the fluid phasē whereC i is the phase average of the concentration of species i in the pore space and q i the charge (in C) of the ionic species i (q i = 0 for neutral species) (see Pride 1994 for the definition of the phase average). This approximation does not mean that we consider that the double layer is thick, but only that the charge distribution is averaged over the pore space. Using this approximation, Revil & Linde (2006) and Revil (2007) have presented recently a model to study electrokinetic properties of deformable porous materials saturated by a viscous (Newtonian) solvent like an aqueous electrolyte. This model was validated through experiments at different scales (Revil et al. 2003;Bolève et al. 2007;Crespy et al. 2008) and generalized to unsaturated conditions by Linde et al. (2007) and Revil et al. (2007). It has been successfully used recently inside a stochastic framework to determine the flow pattern in geothermal fields (Jardani & Revil 2009). The charge densityQ V corresponds to the effective charge density of the diffuse layer per unit pore volume that can be carried by the flow of the pore water relative to the mineral framework. The surface charge densityQ V V p /S (V p is a pore volume and S a surface area of the solid-fluid interface) counterbalances the surface charge density of the mineral surface plus the surface charge density of the Stern layer of sorbed counterions (Leroy et al. 2007). The mean charge density of the diffuse layerQ V can also be related to the total volumetric charge density Q V (traditionally determined by cation exchange capacity measurements in water saturated rocks, see Revil et al. 1998) bȳ where f Q is the fraction of charges located in the Stern layer. The Stern layer is a layer of sorbed counterions on the mineral surface. The surface charge density of the mineral surface is counterbalanced by the surface charge density of the stern plus diffuse layers. The counterions are charges of opposite site to the charge of the mineral surface while the coions have charge of the same charge than the mineral surface. According to the electrochemical models developed by Leroy & Revil (2004) and Leroy et al. (2007) for clay minerals, approximately 90 per cent of the counterions are usually located in the Stern layer in the case of Illite and smectite and over 98 per cent in the case of kaolinite. For silica, an electrochemical model has been discussed recently by Leroy et al. (2008). This model shows that the fraction f Q of counterions located in the Stern layer at the surface of silica grains depends strongly on the salinity (and pH) of the solution. At low salinities (≤10 −3 mol L −1 ), most of the counterions are located in the diffuse layer while at high salinity (>10 −3 mol L −1 ) the counterions are mainly located in the Stern layer. The situation is different for clays for which f Q is much less sensitive to the salinity of the free pore water. Fig. 3 shows that for a wide range of porous media and for electrolytes of different ionic strengths, the excess of charge per unit pore volumeQ V can be related to the DC-permeability k 0 . While most of the data reported in Fig. 3 concern aqueous electrolytes, a similar trend is expected for wet oils as discussed further at the end of this section.
The boundary value problem for the fluid flow is expressed on an averaging disc (see Pride 1994) by where S represents the solid-pore fluid interface andẑ is the unit vector normal to the disc face of the representative elementary volume. This boundary-value problem is exactly the same as the boundary-value problem studied by Pride (1994) for a poroelastic porous body saturated by a Newtonian fluid except that η f should be replaced by η * f . We can therefore easily generalize the results obtained by Pride (1994) to a poroelastic material saturated by a viscoelastic solvent. This yields the following modified Darcy equation is the filtration displacement of the fluid phase relative to the mineral framework (Darcy velocity), φ is the interconnected porosity,ẇ(ω) = φv(ω) is the Darcy velocity, k(ω) is the dynamic permeability, c(ω) is the dynamic coupling coefficient described below and σ is the electrical conductivity of the porous material. The electrical conductivity is given by upscaling the local Nernst-Planck equation. We assume that the wet oil carries a net charge per unit volume to compensate the fixed charge density at the surface of the minerals. In this case, using the Donnan model described by Revil & Linde (2006), the total current density j (in A m −2 ) of the porous material is given by j = σ E +Q Vẇ . In this equation, the last term corresponds to the source current density given by the product between an effective charge density and the Darcy-velocity. This equation differs from the more conventional form, which employs the product of a cross-coupling coefficient with the gradient of fluid pressure and an inertial coupling term. This approach has much less parameters than the classical approach because the charge densityQ V can be determined from the DCpermeability alone. At low frequency, the electrical field E is related to the electrical potential ψ by E = −∇ψ to satisfy ∇ × E = 0. The low-frequency coupling coefficient is given by This equation has been validated on core samples saturated by electrolytic solvents (Revil et al. 2005). The electrical conductivity of the porous material is given by (Johnson 1986) (2007) and Jardani et al. (2007). For porous materials filled with wet-oil, we expect that a similar trend occurs but the volumetric charge density in the pore space would be smaller by comparison with an electrolyte-filled porous medium having the same permeability.
where F is the electrical formation factor defined by the ratio between the tortuosity and the porosity, S is the specific surface conductivity associated with the electromigration of charge carriers (counterions and protons) in the electrical double layer,σ f is the effective conductivity of the pore fluid including the effect of the electrical diffuse layer as discussed by Revil & Leroy (2004) and Revil et al. (2005) and is the characteristic pore size defined by Johnson (1986). A model of S including the contribution of the Stern layer has been developed recently by Leroy et al. (2008) for aqueous solutions and can be used to assess the frequency dependence of this parameter. For simplicity, however, we assume that S , and hence conductivity, are independent of frequency in the developments that follow.
The generalized Darcy's and Ohm's laws appear therefore as cross-coupled constitutive equations whereF f is the applied body-force acting on the fluid phase.
The normalized dynamic permeability is given bỹ where k 0 = 2 /(8F) (Johnson 1986) is the permeability at low frequencies ω Min(ω c , ω m ), ω m = 1/τ m and ω c = η f /(k 0 ρ f F). Both F and are two textural parameters defined by (Johnson 1986) where dV p denotes an integration over the pore space and dS denotes an integration over the pore-solid interface. The normalized electrical field e = ∇ is obtained by solving the following canonical Laplace problem for an average disc in the direction of the macroscopic electrical field (see Pride 1994): has units of length, satisfies Laplace equation ∇ 2 = 0 throughout the pore space with boundary conditions = 0 on z = 0 and = L on z = L, andn · ∇ = 0 on both the internal interface between the solid and fluid phases (n is the unit vector at the solid-fluid interface pointing inside the solid) and on the external surface of the averaging disc (n is the outward unit vector). This approach assumed implicitly the continuity of the solid phase at the scale of a representative elementary volume of the porous rock.
The dynamic streaming potential coupling coefficient is given by The normalized coupling coefficient is given bỹ Therefore, in our approach, the normalized dynamic coupling coefficient is equal to the normalized dynamic permeability, the normalization being made with respect to the DC value of these parameters.
Using the thin double layer assumption and the Boltzmann distribution for the concentrations in the diffuse layer (small size of the diffuse layer with respect to the size of the pores), the dynamic relative coupling coefficient of a porous material saturated by a viscous Newtonian fluid is given by (Pride 1994) and where ζ is the so-called zeta potential (the inner potential of the electrical diffuse layer).
Using the dynamic viscosity defined by eq. (16) and using eq. (47), the dynamic relative coupling coefficient of a porous material saturated by a generalized Maxwell fluid in the thin double layer assumption is given bỹ In eq. (48), the effect of the viscosity appears twice: in the numerator because of the dependence of c 0 with the viscosity and in the last term in the denominator because of the dependence of the critical frequency ω c with the viscosity. It is interesting to see that eq. (48) has the same asymptotic limit at low frequencies as eqs (42) and (46) using the Donnan assumption. This limit is given by 1 + i(ω/ω c ) for a Newtonian fluid.
The excitation of the pore fluid in the resonance frequency band of the viscoelastic fluid yields an amplification of the value of the streaming potential coupling coefficient by several orders of magnitude with respect to the value of the streaming potential coupling coefficient at low frequencies. We expect that the determination of c(ω) in a wide range of frequencies (e.g. 0.01-10 kHz in the laboratory) could help to characterize the properties of the oil contained in the pores through the determination of both τ m , the peak of the relaxation time, and the exponent c, which characterizes the broadness of the distribution of the relaxation times associated with the composition of the oil in terms of polymers.
We can determine also an order of magnitude of the DC-value of the streaming potential coupling coefficient of a sandstone filled with a polar oil. For the Berea sandstone (permeability ∼10 −13 m 2 , porosity ∼0.18), an analysis of the data displayed by Alkafeef et al. (2006) for the streaming current density versus the fluid pressure gradient yieldsQ V = 10 −4 C m −3 . Alkafeef et al. (2006) used a crude oil containing asphaltene, a strongly polar polymer responsible for a high dielectric constant for this fluid. The previous analysis implies a charge density 6 orders of magnitude smaller than existing in a sandstone of same permeability filled by an electrolyte (Fig. 3).
The shear viscosity of the crude oil is equal to 5 × 10 −3 Pa s at 25 • C (Alkafeef et al. (2006). Using σ = 2 × 10 −10 S m −1 , we find a streaming potential coupling coefficient on the order of c 0 = −4 × 10 −5 V Pa −1 . So the low value of the charge density is more than compensated by the effect of the high value of the electrical resistivity of the fluid. Indeed, −4 × 10 −5 V Pa −1 is a very high value by comparison with the value of the streaming potential coupling coefficient of a typical soil or reservoir rock saturated by a brine (typically in the range from -(10 −6 to 10 −8 ) V Pa −1 , see Revil et al. 2003). This implies than even much below the relaxation frequency of the resonance of the solvent, the streaming potential coupling coefficient of a sandstone filled with such an oil can reach a very large value.

T H E M E C H A N I C A L E Q UAT I O N S
The bulk stress tensor is defined bȳ whereT s andT f are the phase average stress tensors in the solid phase and viscoelastic fluid, respectively. The confining pressure and the deviatoric stress of the porous body are defined by wherep s is the mean pressure in the solid phase. The deviatoric stress is then given bȳ whereT D s =T s +p s I is the mean deviatoric stress of the solid phase andπ is the mean deviatoric stress in the pore fluid phase (see Appendix A).
Performing a volume-average of the force balance equations in each phase, eqs (17) and (25), yield a total force balance equation in the time domain where the bulk density and the bulk force are defined by respectively. An expression of the force in terms of source mechanism will be given in Section 6.4. The total stress tensor and the pore fluid pressure are given by (see Appendix A) where the circled cross stands for the Stieltjes convolution product, K U and G U are the undrained bulk and shear moduli of the porous medium, and C and M are two Biot moduli.

T H E M A X W E L L E Q UAT I O N S
As shown by Pride (1994), the local Maxwell equations can be volume-averaged to obtain the macroscopic Maxwell equations (see also Revil & Linde 2006, who used the Donnan model discussed in Section 3). With the Donnan model, the Maxwell equations are where H is the magnetic field, B is the magnetic induction and D is the displacement vector. These equations are completed by two electromagnetic constitutive equations: D = εE and B = μH where ε is the permittivity of the medium and μ is the magnetic permeability. If the porous material does not contain magnetized grains, these two material properties are given by (Pride 1994) where F is the electrical formation factor defined by eq. (43), ε f and ε s are the dielectric constants of the pore fluid and the solid, respectively, and μ 0 is the magnetic permeability of free space. Note that only two textural properties and F are required to describe the influence of the topology of the pore network upon the material properties entering the transport and electromagnetic constitutive equations.
The coupling between the mechanical and the Maxwell equations occurs in the current density, which can be written in the time domain as where j S is the source current density of electrokinetic nature. The body forceF f can be responsible for electromagnetic signals that are directly associated with the source.

Field equations
In this section, we are interested to solve the seismoelectric problem for poroelastic media saturated by Newtonian or generalized viscoelastic fluids. As the case of a poroelastic body saturated by a Newtonian fluid is a special case of the theory develop above, we develop the general field equations below and we show that these equations yields the correct estimates in the limit of a Newtonian fluid. For the seismoelectric problem and neglecting the electrosmotic contribution in the Darcy velocity, the three equations to solve are given, in the time domain, by Darcy's law plus combined equations of motion for the solid and the fluid phases, where M is one of the Biot' moduli (see Appendix A). eq. (67) is derived from eq. (41) in Appendix B, b(t) = η * f (t)/k 0 is defined in Appendix C, η * f (t) is the inverse Fourier transform of η * f (ω) defined by eq. (16),ρ f is defined by eq. (B3). Eq. (68) follows from eq. (54) while eq. (69) follows from eq. (A41) using the relationship α = C/M.
Taking the time derivative of eq. (69) and using Darcy's law in which we have neglected the electroosmotic contribution (the last term of eq. 38) because we are only interested in this paper by the seimoelectric coupling, we obtain, in the frequency domain, the following hydraulic diffusion equation for the pore water (see Appendix B) The assumption that the electroosmotic component can be neglected if we are dealing with the electromagnetic field associated with hydromechanical disturbances has been investigated by Revil et al. (1999) among others. The first term of the left-hand-side of eq. (70) corresponds to the storage term while the second term corresponds to the divergence of the Darcy velocity. The term on the right-hand side of eq. (70) corresponds to a source term for this PDE.
For a 2-D problem, eqs (67) and (68) can be combined to give (see Appendix B) wherep f is related directly to the components of the displacement of the solid and the components ofw by eq. (69). This form is suitable for the implementation into a finite element code. We use Comsol Multiphysics 3.4 adding perfect matched layer (PML) boundary conditions to the problem for the poroelastic waves and associated electrical signals (see Zeng et al. 2001;Jardani et al. 2010). This is to avoid spurious reflections at the side and bottom boundaries of the system. eq. (71) is solved first for the poroelastic propagation problem, which gives the solution for the unknown variables (u x , u z , w x , w z ). The electromagnetic problem is solved in the quasi-static limit. To simplify the problem, we consider that the reservoir is close enough from the sensors (antennas, non-polarizing electrodes and/or magnetometers). In this case, we can neglect the time required by the electromagnetic disturbances to diffuse between the reservoir and the electromagnetic sensors. With this additional assumption, we can model the problem by solving only the quasi-static electromagnetic problem (rather than the low-frequency diffusive problem) for the electrical potential and the magnetic fields where ψ is the electrostatic potential (E = −∇ψ in the quasi-static limit of the Maxwell equations), and j S =Q Vẇ .

Dilatational wave speeds
Combining eqs (69) and (70) yields These two equations have the same structure as eqs (183) and (184) of Morency & Tromp (2008). As shown in Appendix C, the frequency or time dependence of the material properties is however quite different. While Morency & Tromp (2008) investigated dissipation mechanisms associated with the behaviour of the skeleton in terms of an equivalent viscoelastic behaviour of the material, the dissipation mechanism occurring in our model is stemming from the viscoelastic properties of the pore fluid.
In the time domain, the constitutive equations for the total stress tensor and the fluid phase average stress tensor are given in Appendix Ā where the deviatorsd s andd w are defined in Appendix A andd s is the mean strain deviator of the solid phase, and the coefficients G U (t), G U (ω), C G (t), C G (ω), M G (t) and M G (ω) are defined in Appendix A. We use the classical decomposition of the displacements into dilatational components The dilatational wave propagation is obtained by applying the divergence operator to eqs (76) and (77) substituting eqs (78) and (79). This yields where H is the stiffness coefficient of Biot theory (H = K U + (4/3)G U ). We assume plane wave propagation in the x-direction where ω is the angular frequency, l is the complex wavenumber, and z p = ω/l is the complex speed. Using eqs (87) and (88) into (85) and (86), we obtain Eliminating A 1 and A 2 between these two equations yields the following equation for the speed: The two roots of this equation are with z I p > z I I p . These two waves correspond to the complex wave speeds associated with the fast P wave when the solid and the fluid move in phase and the slow dilatational (slow) wave when the solid and the fluid move out of phase. The phase speeds are given by 1/c I p = Re[(z I p ) −1 ] and 1/c I I p = Re[(z I I p ) −1 ] and the inverse quality factors are 1/Q I p = Im(z I p ) 2 /Re(z I p ) 2 and 1/Q I I p = Im(z I I p ) 2 /Re(z I I p ) 2 .

Rotational wave speeds
Using the decomposition of the displacements into rotational components applying the curl operator to eqs (76) and (77), and using eqs (78) and (79), we obtain We assume plane wave propagation in the x-direction where ω is the angular frequency, l is the complex wavenumber, and z S = ω/l is the complex speed. Using these equations, we obtain Eliminating B 1 and B 2 between these two equations yields the following equation for the speed: a z 4 The two roots of this equation are with z I S > z I I S . These two waves correspond to the complex wave speeds associated with a fast rotational S wave and a slow rotational S wave. The slow rotational wave represents the classical shear wave while the fast rotational wave is related to the fact that the fluid can sustain shear stresses above its resonance frequency. This wave does not exist at low frequencies (ω ω m ) because the pore fluid behaves as a Newtonian fluid. The phase speeds of these two waves are given by 1/c I S = Re[(z I S ) −1 ] and 1/c I I S = Re[(z I I S ) −1 ] and the inverse quality factors are 1/Q I S = Im(z I S ) 2 /Re(z I S ) 2 and 1/Q I I S = Im(z I I S ) 2 /Re(z I I S ) 2 , respectively. To check the consistency of the model, we can look for the case where G * f (ω) = 0 and G U = G f r , C G (ω) = 0, and M G (ω) = 0. In this case, we obtain and therefore In absence of dissipation, this second wave corresponds to the classical shear wave with phase speed (see Morency & Tromp 2008).

Mechanical source
The macroscopic source term F is linearly partitioned between the solid and the fluid phases, see eq. (56). It can be written in terms of the moment tensor M by (Dahlen & Tromp 1998) where x S refers to the position of the source, S(t) is the source-time function, and δ(x) is the Dirac (delta) function. This source term arises in eq. (71) to determine the initial values of the displacements of the solid phase and the relative displacement of the fluid phase with respect to the solid. Then the relative displacement of the fluid phase with respect to the solid phase creates a source term in the electrostatic and magnetostatic field equations, eqs (74) and (75).

S Y N T H E T I C C A S E S T U D I E S
Experiment #1 corresponds to the simulation of the seismoelectric response of an oil-filled reservoir (see Fig. 4) in an half-space with the oil considered as a Newtonian fluid. The ground surface corresponds to the top surface of the system where the normal component of the electrical field vanishes. The material properties of the embedding medium (with the properties of a sandstone) and the oil reservoir are reported in Table 1. We consider a volumetric source located 50 m below the ground surface generating P waves. The reservoir is 100 m thick and its top surface is located 250 m below the ground surface (Fig. 4a). The time dependence of the source S(t) is a Ricker with a dominant frequency of 30 Hz. We consider the arrival of the seismic and seismoelectric signals at a station P 1 (−50 m, 0 m) corresponding to a single electrode. The reference for the electrical potential recorded at this electrode is located at position Ref(−300 m, 0 m) (see Fig. 4). The four edges are absorbing boundaries for which we use PML boundary conditions (see Jardani et al. 2010 for more details). To solve numerically the problem, we use the finite element modelling Comsol Multiphysics 3.5 to simulate both seismograms and electrograms at the ground surface.
Both seismoelectric conversion and coseismic electrical signals are generated at the reservoir boundaries (Figs 4b, 5 and 6). Fig. 6 shows two electrograms for station P1. The first type of signals corresponds to the co-seismic electrical signal associated with the propagation of the P wave. It is labelled 'Coseismic' in Figs 5(a) and 6. Other coseismic signals are associated with the reflected P waves at the top and bottom interfaces of the reservoir. These coseismic signals are labelled RCS1 and RCS2, respectively. These coseismic There is a signal associated with the source itself that diffuses nearly instantaneously at the observation station. We decided not to model this signal. In the early times of the electrogram, various contributions include the seismoelectric conversions of the P wave reaching the top of the reservoir (2) (labelled IR1) and the bottom of the reservoir (3) (labelled IR2). In addition, coseismic signals are associated with the direct wave (1) (labelled 'coseismic') and the reflected waves (4) and (5) (labelled RCS1 and RCS 2). signals occur when a seismic wave travels through a porous material, creating a relative displacement between the pore water and the solid phase. The associated current density is balanced by a conduction current density. It results an electrical field travelling at the same speed as the seismic wave. The second type of seismoelectric signals correspond to converted seimoelectric signals associated with the arrival of the P waves at the top and bottom interface of the oil reservoir. These converted seimoelectric signals are labelled IR1 and IR2, respectively. When crossing an interface between two domains characterized by different properties, a seismic wave generates a time-varying charge separation, which acts as a dipole radiating electromagnetic energy. These dipoles oscillate with the waveform of the seismic waves. Because the electromagnetic diffusion of the electrical disturbance is very fast (instantaneous in our simulations), the seismoelectric conversions are observed nearly at the same time by all the electrodes but with different amplitudes. The seismoelectric conversions appear therefore as flat lines in the electrograms shown in Fig. 5(a) Experiment #2 uses the same geometry as Experiment #1 but the fluid in the reservoir is now a viscoelastic wet-oil. The dynamic permeability response of the reservoir filled with the viscoelastic fluid is shown in Fig. 7(a). The time function of the source is a Dirac and we investigate the electrical field response in the frequency band (1-1000 Hz). A plot of the vertical component of the electrical field at an observation point at depth is shown in Fig. 7(b). The maximum at the electrical field occurs at the same frequency that the resonance of the viscoelastic fluid. The amplification of the electrical field reaches several orders of magnitude with respect to the DC-value. This very strong amplification of the signal could be the basis of a new detection and mapping method of heavy oils and DNAPL in the ground.

C O N C L U D I N G S TAT E M E N T S
We have developed a model to simulate the seismoelectric coupling in a poroelastic body saturated by a viscoelastic or a Newtonian fluid ignoring the electroosmotic term in the Darcy equation and assuming that the electromagnetic disturbances generated by the seimoelectric conversions move instanteanously though the conductive Earth. The following conclusions are reached.
(1) When the pore fluid is a wet oil, it can be also considered as a viscoelastic non-aqueous solvent with a finite electrical conductivity and a relatively high dielectric constant (but always lower than the dielectric constant of water).  Fig. 4). The strongest signal on the electrogram corresponds to the coseismic disturbance associated with the direct wave (see Fig. 4). RCS1 and RCS2 stand for the coseismic disturbances associated with the reflected P waves (see Fig. 1). IR1 and IR2 stand for the two seismoelectric disturbances associated with the seismoelectric conversions at the top and bottom of the reservoir.   Fig. 4. The fluid of the reservoir is viscoelastic. The source function is a Dirac and we investigate the electrical response in the frequency band (1-100 Hz). The dynamic permeability versus frequency shows the relaxation peak associated with the resonance of the viscoelastic fluid. The peak of the vertical component of the electrical field E at a remote self-potential station located at x = 150 m and y = 0 m.
(2) Resonance of this solvent occurs in a spectrum of frequencies that depends on its composition. The distribution of the relaxation times of such a fluid can be described by a Cole-Cole distribution. For heavy oils and DNAPL, the resonance frequencies belongs to the seismic band used by conventional seismic methods (10-100 Hz).
(3) If the porous body is excited by seismic waves with frequencies inside this frequency band of resonance of the viscoelastic fluid, the current density associated with the movement of the viscoelastic solvent can be very strong with an enhancement of the streaming potential coupling coefficient of several orders of magnitude with respect to its value at low frequencies. In all the spectrum of frequencies, the source current density, in turn, creates electromagnetic disturbances that can diffuse away from the source of current density. If we are close enough to the source volume of current, the electromagnetic problem can be considered in the quasi-static limit of the Maxwell equations.
(4) A new type of shear wave has been discovered. This new shear wave is controlled by the mechanical properties of the viscoelastic fluid and is related to the resonance of the fluid and its ability to bear shear stresses. This resonance creates a resonance in the electrical field of electrokinetic nature that may be strongly amplified by several orders of magnitudes. Future works will be dedicated (1) to the check of the present laboratory in the laboratory and in the field and (2) to the extension of the theory to the case where the pore space of a porous material is filled with two immiscible fluid phases like water and oil looking at the effect of the wettability of the oil phase upon the seismoelectric response.

A C K N O W L E D G M E N T S
We thanks Terry Young for his support at CSM. We thank Kasper van Wijk and one anonymous referee for their very constructive reviews and Mark Everett for his work as Associate Editor.