Seismoelectric effect in Lamb’s problem

Seismoelectric effect is studied in the framework of classical Lamb’s problem with impulse or timevariable mechanical action on an elastic porous half-space. Radiation of elastic waves gives rise to pressure variation of groundwater fluid contained in pores and cracks. This causes the generation of telluric electric fields and currents due to the seismoelectric effect. A diffusion type equation is applied to describe the variations of the pore pressure and telluric electric field. Particular emphasis has been placed on the properties of seismoelectric signals caused by Rayleigh wave propagation since this wave has maximal amplitude at a considerable distance from the seismic source. For practical purposes and geophysical application, the co-seismic phenomena related to seismoelectric effect are examined in more detail.


Introduction
In the field of applied seismology the Lamb's problem [Lamb, 1931] is of great utility especially in the modelling of seismic wave originated from earthquakes, explosions, meteorite impacts on the Earth and etc. In the Lamb's problem, the source of acoustic waves is generally assumed to be an impulse mechanical action on the elastic halfspace. There are several types of such waves including longitudinal, transverse and surface waves, which can propagate in elastic media at different speeds. Actual geophysical media consist of a great variety of rocks, which contains inhomogeneous inclusions such as quartz, calcite and dolomite as well as cracks, pores filled with groundwater and etc. The seismic wave propagation in the ground is accompanied by the deformation of the pore space which results in the pore vibrations, diffusion of underground pore fluid followed by the seismoelectric phenomena and weak perturbations of the geomagnetic field [e.g., Surkov and Hayakawa, 2014].
The co-seismic electromagnetic variations are usually observed during the passage of seismic waves through the ground-based magnetometers and antennas or buried electrodes [Ivanov, 1940;Martner and Sparks, 1959;Eleman, 1965;Anisimov et al., 1985;Iyemori et al., 1996]. There are two basic physical mechanisms of these coseismic phenomena: the seismoelectric effect in porous water-saturated media and the perturbation of the Earth's magnetic field due to the motion of the conductive ground [Surkov et al, 2018]. The seismoelectric effect builds up as a result of electric charge separation caused by the electro-chemical processes in the electrolyte solution, which is contained in the underground fluid [Frenkel, 1944]. The surfaces of cracks and pores can adsorb ions of certain polarity (more frequently the negative-charged ions) from the electrolyte solution whereas the underground fluid accumulates the excess ions of opposite/positive sign. The seismic wave induces the underground fluid motion, which results in the generation of electrokinetic currents due to entrainment of the positive ions with the fluid. The electrokinetic current closes through a conduction current flowing in the fluid and conductive rocks. Notice that the seismoelectric effect manifests itself in electric field predominantly rather than in magnetic field.
The co-seismic variations of the telluric electric potentials have been observed after the occurrence of earthquakes with magnitude ≥ 5 at short epicentral distances Mogi et al., 2000]. Analysis of the experimental data has shown that there are two types of co-seismic electric variations measured with the electrodes, which were embedded in the ground. The oscillating perturbations arising at the time of seismic wave arrival can be interpreted in terms of the seismoelectric effect. The other type is the co-seismic variations gradually declining during 1 minute after earthquake occurrence. This attenuation effect has been assumed to be due to the changes in hydrologic status of the ground surface layer . Electric and magnetic perturbations have been observed after the Izmit earthquake ( = 7.6) in Turkey on August 1999 immediately upon arrival of the seismic P-wave at the magnetotelluric stations [Honkura et al, 2002].
Acoustic and electrokinetic phenomena in porous water-saturated media are related through a set of coupled equations for acoustic and electrodynamic parameters [Frenkel, 1944;Pride, 1994;Pride and Garambois, 2005;Hu and Gao, 2011]. When considering the seismic problem, the electroosmotic phenomenon can be neglected. In this case, the basic set of equations can be split into two different sets in such a way that one of these sets includes only the acoustic field variables, independently of the electric field variables. This means that the acoustic variables contained in the second set play the role of known sources of seismoelectric variations. The seismoelectric effect has been studied in more detail for the case of a harmonic acoustic wave propagating in an infinite homogeneous space [e.g., Surkov et al., 2018]. More complex models of the seismic source are necessary in order to study the seismoelectric phenomenon associated with earthquake and explosion. The models need to take into account the impulse character of the seismic source and the influence of boundary conditions at the ground surface. In this study the seismoelectric effect is analyzed on the basis of solution of the classic Lamb's problem; that is, the problem of an impulse or transient impact on an elastic porous half-space.

Perturbation of pore fluid pressure caused by the propagation of the acoustic wave
The electrokinetic effect in porous water-saturated medium arises under influence of the pore fluid pressure gradient ∇ . The field of extrinsic force related to this effect is described by the vector E = -∇ , where is the streaming potential coupling coefficient. The total current density in the porous medium consists of the electrokinetic current density E and the conduction current -∇ , where is the mean conductivity of the porous medium and is the electric field potential. Considering the steady-state distribution of the electric current in porous medium, the potential is described by the following equation We shall restrict our consideration to a homogeneous medium with constant and . Thus equation (1) reduces to The electric potential as well as the normal component of the total current density have to be zero at the boundary between the ground and the atmosphere. In addition, the pore fluid pressure has to be equal to atmospheric pressure on the ground surface. This requirement is replaced by the approximate boundary condition = 0 since the atmospheric pressure is small. In such a case the solution of equation (2) has a simple view (e.g., [Surkov and Hayakawa, 2014]):

Vadim V. Surkov et al.
The equations related the volumetric strain of porous medium and the variation of pore fluid overpressure (with respect to hydrostatic pressure) have been derived by [Frenkel, 1944]. These equations include the second order time-derivatives which can be neglected in the low frequency region. In this approximation the equations are simplified and reduced to the diffusion-type equations: where denotes the coefficient of diffusion of the fluid pressure perturbation in porous medium, stands for the viscosity of the pore fluid while and are the medium porosity and permeability, respectively. Here we made use of the following abbreviations: where is the bulk modulus of compressibility of pore fluid, is the bulk modulus of solid which forms a medium skeleton and is the bulk modulus of dry porous rock (the medium skeleton without the pore fluid).
The porous medium is treated as a homogeneous elastic half-space bounded from above by the atmosphere.
Since the medium porosity is assumed to be small, the velocities of elastic waves and other medium parameters depend only slightly on the porosity. Additionally, it is conjectured that the wave damping due to viscous friction in the fluid can be neglected. The volumetric strain in equation (4) is treated as a given function. In order to find this function we first consider the Lamb's problem for the elastic half-space. Suppose the origin of coordinate system is located on the surface of the half-space while the z-axis is pointed vertically downward. The boundary surface is subjected to the normal force distributed on the surface. The normal component of the stress tensor corresponding to this force is assumed to be a given function at = 0: where ₀ and are the force amplitude and typical time of the force action, ₀ is the characteristic radius of the area subjected to the normal force, = ( ²+ ²)¹ / ² is polar radius and ( ) denotes Dirac delta-function of time. Notice that if ₀ → 0 then ₀ / ( ²+ ₀² )³ / ² → ( )/ . In this extreme case we come to the known axially symmetric Lamb's problem on point normal force acting on elastic half-space.
The cylindrical coordinates z, r are convenient to use since the normal force (6) is axially symmetrically distributed on the surface. The solution of the corresponding boundary-value problem for the elastic half-space enables one to determine the radial and vertical components of the displacements of the elastic medium [Petrashen et al., 1950;Brekhovskikh, 1980;Nikitin et al., 2013]: Here the following designations are used for Laplace integrals where ( ) = (2+ ²)² -4 ₁ ₂ denotes Rayleigh function. Here we made use of the following abbreviations (11) where and are the velocities of longitudinal and transverse elastic waves while is the undisturbed medium density.
The volumetric strain is given by [Landau and Lifshitz, 1970]: Substituting equations (7) and (8) for and into equation (12), yields (13) where (14) Taking into account the construction of the expressions for functions (13) and (14) we seek for the solution of equation (4) in the form: Substituting equations (13) and (15) for and into equation (4) we obtain the following equation for the unknown function : where (17) Here the primes denote derivatives with respect to z. Assuming the pore overpressure is zero at the surface z = 0 and taking into account that the solution of equation (16) has to be limited with z → we obtain Substituting equations (17) and (18) for and into equation (15), we find the overpressure of the underground fluid in pores: As is seen from equation (3) the distribution of the electric potential in the ground can be expressed through via = -. Thus, we have derived the general solution of the above problem for the case of a pulsed seismic source. The solution of the more general problem for an arbitrary source can be derived by convolution of the above solution with the source function. We shall concern this problem in the fourth section.
The contribution of the longitudinal and Rayleigh waves to the variation of the pore fluid pressure is associated with the singularities of the integrand in equation (19) that are discussed in Appendix. The function ₁ associated with the transverse wave does not enter equation (19). The reason is that the transverse wave does not create the volume deformation of the medium and thus it has not effect on both and the seismoelectric signal. However, in the high frequency range the shear horizontal (SH) wave can serve as a source of the seismoelectric effect even if = 0 (e.g., [Han and Wang, 2001;Monachesi et al., 2018;Gao et al., 2019]). In the theory, the electric signal produced by the SH-wave is proportional to a relative acceleration of the fluid phase with respect to the solid one [Han and Wang, 2001]. This source does not enter equation (4) for since this relative acceleration is proportional to the frequency squared and thus is small in the seismic frequency range (below or of the order of several Hz). In the remainder of this paper we shall return to this question.
It should be noted that the derived solutions are valid for the weakly porous homogeneous media where the dissipation and dispersion of acoustic waves can be neglected. In the more general case of dissipative porous watersaturated media, it is necessary to solve a set of coupled equations that includes both the acoustic and electrical parameters [Frenkel, 1944;Pride, 1994].

Variation of pore fluid pressure caused by the Rayleigh wave
In this section we seek for the variation of the pore fluid pressure due to the propagation of the Rayleigh wave.
The contribution of this wave to is defined by the poles of the integrand in equation (19). The coordinates of these poles are = ± . Here = / where denotes Rayleigh wave velocity [Lamb, 1931]. Calculations of the corresponding pole residues are found in Appendix. The result of these calculations is given by Here we made use the following abbreviations The diffusion coefficient of the pore fluid pressure is directly proportional to the medium porosity which can vary within a very wide range. For example, the average porosity of the typical rocks are of the order of (0.2-6)·10⁻³, 0.04 -0.3 and 0.2 -0.3 for granite, gneiss and tuff, respectively [Mavko et al., 2009]. In order to estimate the diffusion coefficient and other parameters we use the following numerical values: = 0.1, = 2 GPa, = 10 GPa and = 20 GPa whence it follows that = 0.714. The water viscosity coefficient varies within (0.6 -1.8)·10⁻³ Pa·s in the temperature interval from 0° C to 40° C. The permeability coefficient of the actual rocks also varies in a very wide range from 10⁻¹⁴ to 10⁻⁹ m 2 [Teodorovich, 1958]. For the permeable and well permeable rocks at =10⁻¹² -10⁻¹⁰ m 2 the diffusion coefficient in equation (4)  Here ₀ is the characteristic scale of the function decrease. In order to estimate this parameter, we use = 1.4 km/s and the above range of values of the diffusion coefficient. Then we get the following estimate: ₀ ≪ ⁻¹ under the requirement that ≫ * = 2 / = 1.2·10⁻³ -0.34 m. In this approximation the first exponential function in the square brackets in equation (20) falls off more rapidly with than does the second exponential function. Therefore, with these parameters the function exp(-) in equation (20) becomes insignificant in this range of depths . In addition, one can neglect the summand in the denominator of equation (20) under the requirement that ≫ / ²) ≈ 2.6·10⁻⁴ -0.08 m. In practice, to measure telluric potentials the electrodes are buried in the ground for a few tens of cm so that both of the above conditions hold true.
This approximation can be taken into account by setting formally = 0 in equation (20). As a result we obtain: where = ₂ + ₀. Performing integration in equation (23) yields [Grandshteyn and Ryzhik, 2007] (24) where ± = ² + ± ² ¹ / ² . (24) near the wave front and far away from the seismic source. By making assumption that | | ≪ and ≪ , equation (24) can be simplified (25) where (26) Let ₁ = / be a new dimensionless variable. Taking derivative of this expression with respect to we get (27) Now let us make some simple estimates to gain a better understanding of the application area of the results obtained above. For this purpose we consider an impulsive pressure source situated at some point in elastic pore space. After the source is switched on, the spherically symmetric P-wave and the diffusion of the pore fluid begin to propagate simultaneously from this point. The front radius of the seismic P-wave increases with time as ~ whereas the size of area covered by the diffusion increases as ~ 2( )¹ / ² (e.g., [Surkov and Hayakawa, 2014]).

Let us introduce a new variable =in order to analyze equation
From these relations it follows that at first the diffusion front propagates faster, that is, ahead of the seismic wave.
The diffusion velocity decreases with time, so the P-wave catches up with the diffusion front at some distance * .
Equating and , one can find this distance: * ~ 2 * ≈ 4 / ≈ 0.2 -70 cm. This value is much smaller as compared to the seismic wavelength. This means that the process of the pore fluid diffusion is insignificant because the pore fluid overpressure is mainly controlled by the variations of the seismic wave parameters.

Seismoelectric effect caused by the Raleigh wave
In order to estimate the Rayleigh wave contribution to the seismoelectric effect, equation (25) for pore fluid overpressure needs to be substituted into equation (3) for telluric potential: The strength of the telluric electric field due to Rayleigh wave is derivable from a potential through = -∇ .
The most interesting is the horizontal component of the telluric field which dominates near the ground surface (29) where (30) The dimensionless dependences of the potential and the radial component of the electric field in a reference frame moving at Rayleigh wave velocity are shown in Figures 1a and 1b as functions of the dimensionless distance ₁. As is seen from these Figures, the potential and electric field reach peak values of the order of and , respectively, in the region around the point ~ . For the surface layer of the ground; that is, under requirement that 0 ≫ 2, one can assume that ~ 0. Taking into account this approximation and substituting ₀ = ₀² ₀ into equations (28) and (29), we obtain (31) Here 0 is the characteristic radius of the area subjected to the force ₀ while ₀ = ₀ / ( ₀² ) is the pressure caused by this force. To take one example, consider the electrokinetic effect resulted from the point explosion on the ground surface.
The approximate dependence of the parameter = / on Poisson coefficient (a coefficient of transverse strain) is given by [Landau and Lifshitz, 1970]:  [Sorochan, 1987].
The dependence of = / on is given by (33) By setting the Poisson coefficient = 0.3 and substituting this value into equations (21), (32) and (33)  other parameters of the medium in which the measuring electrodes are located [Jouniaux et al., 2000;Vinogradov et al., 2010].
The explosion parameters are chosen to be ₀ = 0.1 GPa, = 0.1 s and 0 = 2 m. Additionally we suppose that = 10⁻⁶ -10⁻⁸ V/Pa and ² ≈ = 10 GPa. Substituting these and above-mentioned values of the parameters into equation (31), we obtain that = 1 -100 mV/m at the distance of =10³ km from the explosion site. Typically, the natural electric noise amplitude in the ground is about 1-10 V/m; that is, several magnitudes smaller than .
This means that the seismoelectric signal produced by such kind of explosion can be detected at the background of the telluric noise.

Seismic sources with arbitrary time dependence
To this point the impulsive seismic source has been examined. Now let us suppose that the pressure acting on the ground surface is described by an arbitrary source function ( ). In such a case one must replace time with -, the parameter with ( ) and then make a convolution of the above solution with the function ( ). Hence the function ( , ) in the above solution has to be replaced by the following function (34) where = -( -) = + . By analogy with equations (29) and (30), we can thus write where 1 = / while and / ₁ are given by equations (28)-(30).
By way of illustration let us first consider the pressure pulse described by a step function ( ), which equals unity within the interval 0 ≤ < or zero as > . In such a case it is suitable to perform integration in equation (34). As a result, we obtain ( > ) The co-seismic electromagnetic phenomena have been frequently observed after earthquake occurrence away from the epicenter [e.g., Surkov and Hayakawa, 2014 and references herein]. It is usually the case that the far-field displacement components of the seismic wave have a shape largely similar to the damped oscillations. To illustrate this property of seismic wave, let us now approximate the source function ( ) as follows where is the seismic wave frequency and denotes the characteristic decay time.
Substituting equation (41) for ( ) into equation (35)  From Figure 3, the typical frequency of the electric field vibrations is close to though the electric signal is phase shifted with respect to the seismic variations of the pressure or mass velocity, which are described by the function ( ). Such co-seismic geoelectric field changes after earthquakes occurrence have been observed by many researchers [e.g., Eleman, 1965;Nagao et al., 2000;Honkura et al., 2002;Balasco et al., 2014;Romano et al., 2018].

Discussion and Conclusion
The general solution given by equations (3) and (19) shows seismoelectric effects associated with propagation of P-wave and Rayleigh wave. The transverse wave cannot excite such effects in the low frequency region. The reason is that transverse waves such as the S-wave or Love wave do not create a volume strain and thus they have not influence upon the pressure gradient of pore fluid. However, the observations have shown that the co-seismic electromagnetic signals can be excited by both P-and S-waves [e.g., Balasco et al., 2014;Romano, 2019]. Our rationale is that these co-seismic signals can be interpreted in terms of two main underlying mechanisms: (a) electrokinetic effect in porous medium [e.g., Frenkel, 1944;Martner and Sparks, 1959;Pride, 1994;Hu and Gao, 2011;Surkov and Pilipenko, 2015] and (b) a geomagnetic inductive (GMI) effect [e.g., Knopoff, 1955;Guglielmi, 1986;Gorbachev and Surkov, 1987]. The GMI perturbations during the seismic wave propagation across an observation site arise from the generation of electric currents in the conductive layers of the ground due to the Vadim V. Surkov et al.
10 ground motion in the geomagnetic field. Thus, one may assume that the GMI mechanism can explain the electromagnetic perturbations observed during the S-waves propagation.
The main difference between the seismoelectric and GMI mechanisms is that the GMI disturbance should spread from a seismic wave front along a conductive crust in a diffusive way. Depending on the ground conductivity the velocity of the GMI diffusion can even supersede the seismic wave velocity, thereby producing an "electromagnetic precursor" of a seismic wave front [Surkov 1997]. Such a kind of electromagnetic changes in both magnetic [Honkura et al., 2002] and electric components [Balasco et al., 2014;Romano, 2018] before the arrival of seismic waves has been occasionally observed. However, in practice, it is difficult to conclude unambiguously which of those mechanisms prevails because of a large variability of the realistic crust parameters, such as porosity, permeability, conductivity and etc. [Surkov et al. 2018].
In order to estimate the contribution to the co-seismic signal from the seismoelectric effect, one need information on the distribution of the streaming potential coupling coefficient around an observation point. In this study, we have shown that the coefficient has a much greater influence on the seismoelectric signals as compared with the effect of the fluid diffusion coefficient .
Away from the seismic source, the Rayleigh wave has the greatest amplitude as compared to other seismic waves, and thus it can provide a considerable contribution to the seismoelectric effect. Analysis of the solutions has shown that the electric field perturbations are mainly concentrated near the seismic wave front and propagate in the medium at the speed of Rayleigh wave. In the case of the pulsed source, the seismoelectric signal has a shape somewhat similar to a bipolar pulse unlike the seismic signal. As the seismic source time function resembles a step function, the seismoelectric perturbations are mainly grouped near the front area and in the vicinity of a trailing edge of the seismic signal, i.e. in those areas where the pressure gradient of the pore fluid has the greatest value.
The seismoelectric effect caused by surface detonation has been estimated too. It was shown that for a pressure amplitude in the source of about 0.1 GPa, the electric perturbations in the ground can exceed the background telluric variations at distances up to thousand kilometers.
Modelling of the seismoelectric effect caused by a seismic wave coming from the earthquake focal zone has been performed by using the source function in the form of damped oscillations. The co-seismic electromagnetic signals are shown to be excited practically synchronously with the passage of seismic waves through an observation point though the electric signal is phase shifted with respect to the seismic signal, which is consistent with the observations. The amplitude of the seismoelectric signals has been shown to decrease with distance in the same fashion that the pressure and other Rayleigh wave parameters fall off, i.e. inversely proportional to ¹ / ². However, the exact values of the amplitudes of the seismoelectric signals cannot be predicted in advance because of the uncertainty of the streaming potential coupling coefficient, which can vary widely depending on the porosity and permeability of actual rocks, the mineral composition of groundwater and other parameters of porous media.