Layout 6 1 ANNALS OF GEOPHYSICS, 63, 4, SE440, 2020; doi:104401/ag-8256 Seismoelectric effect in Lamb’s problem Vadim V. Surkov*,1,2, Valery M. Sorokin1 and Аleksey K. Yaschenko1 (1) Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation of the Russian Academy of Sciences, Moscow, Troitsk, Russia (2) Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia Article history: received June 30, 2019; accepted February 17, 2020 Abstract Seismoelectric effect is studied in the framework of classical Lamb’s problem with impulse or time- variable 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. Keywords: Seismology, Ground motion, Waves and wave analysis, Groundwater processes, Fluid geochemistry. 1. 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 half- space. 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 co- seismic 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 [Nagao et al., 2000; 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 [Nagao et al., 2000]. 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. 2. 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 ∇ · 𝜎 (E�� ‒ ∇𝜑) = 0. (1) We shall restrict our consideration to a homogeneous medium with constant 𝜎 and 𝐶��. Thus equation (1) reduces to ∇²(𝜑 + 𝐶��𝑃� ) = 0. (2) 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]): 𝜑 + 𝐶��𝑃� = 0. (3) Vadim V. Surkov et al. 2 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: (4) 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: (5) 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: (6) 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]: (7) (8) Here the following designations are used for Laplace integrals (9) 𝜕𝑃� 𝜕𝑡 = 𝐷𝛻²𝑃� ‒ 𝛼 , 𝐷 = , 𝜕𝜃 𝜕𝑡 𝑛𝛽𝑘�𝐾� 𝜂 𝛼 = (𝐾� ‒ 𝐾)𝛽, 𝛽 = ,𝐾�𝐾� (𝐾� ‒ 𝐾)𝐾� + 𝑛𝐾� (𝐾�‒𝐾�) 𝜎�� = ,𝐹₀𝑇𝛿(𝑡)𝑟₀ 2𝜋(𝑟²+𝑟₀² )³/² ∞ ₀𝑢� = 𝐴� ℎ� �𝑞,𝑧,𝑡�exp�‒𝑞𝑟₀�𝑞J₁�𝑞𝑟�𝑑𝑞, ∞ ₀𝑢� = 𝐴� ℎ� �𝑞,𝑧,𝑡�exp�‒𝑞𝑟₀�𝑞J₀�𝑞𝑟�𝑑𝑞. 𝜀+𝑖∞ 𝜀‒𝑖∞ ℎ� = � ��2 + 𝜁²�exp�‒𝑧𝑞𝑠₂�‒2𝑠₁𝑠₂ exp�‒𝑧𝑞𝑠₁�� ,exp (𝑞𝑡𝐶�𝜁 )𝑑𝜁 𝑅(𝜁) 3 Seismoelectric effect in Lamb’s problem Vadim V. Surkov et al. 4 (10) 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]: (12) 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: (15) Substituting equations (13) and (15) for 𝜃 and 𝑃� into equation (4) we obtain the following equation for the unknown function 𝑝: (16) 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 (18) Substituting equations (17) and (18) for 𝐵 and 𝑝 into equation (15), we find the overpressure of the underground fluid in pores: (19) 𝜕𝑢� 𝜕𝑧𝛳 = �𝑟𝑢�� + . 1 𝜕 𝑟 𝜕𝑟 ∞ ₀𝛳 = 𝐴�ℎ� �𝑞,𝑧,𝑡�exp�‒𝑞𝑟₀�𝑞J₀�𝑞𝑟�𝑑𝑞, 𝜀+𝑖∞ 𝜀‒𝑖∞ ℎ� = � ��2 + 𝜁²�𝜁²exp�‒𝑧𝑞𝑠₂� .exp (𝑞𝑡𝐶�𝜁 )𝑑𝜁 𝑅(𝜁) 𝑞 𝛾² ∞ ₀𝑃� = � � � 𝑝�𝑞,𝑧,𝜁�exp�𝑞𝑡𝐶�𝜁�𝑑𝜁�J₀�𝑞𝑟�𝑑𝑞. 𝜀+𝑖∞ 𝜀‒𝑖∞ 𝑝˝ ‒ 𝜆²𝑝 = ‒ exp�‒𝑧𝑞𝑠₂ ‒𝑞𝑟₀�,𝐵𝑞³𝜁³(2+𝜁²) 𝐷𝑅(𝜁)𝛾² 𝐵 = , 𝜆²= 𝑞² + .𝛼𝐹₀𝑇 4𝜋²𝑖𝜌 𝑞𝜁𝐶� 𝐷 𝜀+𝑖∞ 𝜀‒𝑖∞ ℎ� = � ��2 + 𝜁²�exp�‒𝑞𝑧𝑠₂�‒2exp�‒𝑞𝑧𝑠₁�� ,𝑠₂exp (𝑞𝑡𝐶�𝜁 )𝑑𝜁 𝑅(𝜁) ∞ ₀𝑃� = 𝐵� � � exp�𝑞𝑡𝐶�𝜁‒𝑞𝑟₀�𝑑𝜁� 𝑞²J₀�𝑞𝑟�𝑑𝑞. 𝜀+𝑖∞ 𝜀‒𝑖∞ 𝜁²(2+𝜁²){exp(‒𝜆𝑧)‒exp(‒𝑧𝑞𝑠2)} 𝑅(𝜁)(𝑞𝐷𝜁‒𝐶�𝛾²) 𝐴 = , 𝑠₁ = �1 + 𝜁²�¹/² , 𝑠₂ = �1 + �¹/² , 𝛾 = ,𝐹₀𝑇 4𝜋²𝑖𝜌𝐶� 𝜁² 𝛾² 𝐶� 𝐶� 𝑝 = exp�‒𝑞𝑟₀�.𝐵𝑞²𝜁²(2+𝜁²){exp(‒𝜆𝑧)exp(‒𝑧𝑞𝑠2)} 𝑅(𝜁){𝑞𝐷𝜁‒𝐶�𝛾²} 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 water- saturated media, it is necessary to solve a set of coupled equations that includes both the acoustic and electrical parameters [Frenkel, 1944; Pride, 1994]. 3. 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 (20) Here we made use the following abbreviations (21) 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⁻⁹ m2 [Teodorovich, 1958]. For the permeable and well permeable rocks at 𝑘� =10⁻¹² ‒10⁻¹⁰ m2 the diffusion coefficient in equation (4) is estimated as 𝐷 ≈ (0.008‒2.4)·10⁻² m2/s. Using this numerical estimate one can derive the approximate expression for the integral in equation (20). At first, notice that the real part of the exponent of exponential function exp(‒𝑧𝑞𝜒) is smaller than unit under requirement that (22) 𝑞<𝑞₀ = 𝑧⁻¹�1+ �⁻¹/².𝑧²𝐶�² 4𝐷² 𝑃�,� = �𝑞²exp�‒𝑞𝑟₀�J₀�𝑞𝑟� Re � �exp(‒𝑧𝑞𝜒)‒exp(‒𝑧𝑞𝑝₂)��𝑑𝑞.𝛼𝐹₀𝑇(2 ‒ 𝜉²)𝜉 4𝜋𝜌𝑄(𝜉) exp(𝑖𝑞𝑡𝐶�) 𝑞𝐷𝜉+𝑖𝐶�𝛾² ∞ ₀ 𝑄(𝜉) = 2‒𝜉² ‒ ‒ , 𝜒 = �1+ �¹/², 𝑝₁ = �1 ‒ 𝜉²�¹/², 𝑝₂ = �1 ‒ �¹/².𝜉² 𝛾² 𝑖𝐶� 𝑞𝐷 𝑝₂ 𝑝₁ 𝑝₁ 𝑝₂𝛾² 5 Seismoelectric effect in Lamb’s problem 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: (23) where 𝜇 = 𝑧𝑝₂ + 𝑟₀. Performing integration in equation (23) yields [Grandshteyn and Ryzhik, 2007] (24) where 𝑄± = �𝜇² + �𝐶�𝑡 ± 𝑟�²�¹/². Let us introduce a new variable 𝑥 = 𝑟 ‒ 𝐶�𝑡 in order to analyze equation (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]). 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. 𝑃�,� = ‒ , 𝛼𝐹₀𝑇(2 ‒ 𝜉²)𝐶� 2⁵/²𝜋𝜌𝑄(𝜉)𝐶�² 𝑑² 𝑑𝜇² {𝐶�² 𝑡² ‒ 𝑟² ‒𝜇²+𝑄₊𝑄₋}¹/² 𝑄₊𝑄₋ 𝑃�,� = ‒ 𝑓(𝑥,𝜇) , 𝛼𝐹₀𝑇(2 ‒ 𝜉²)𝐶� 8𝜋𝑟¹/²𝜌𝑄(𝜉)𝐶�² 𝑓(𝑥,𝜇) = .�(𝜇² + 𝑥²)¹/² ‒𝑥�¹ /² (𝜇²+𝑥²)¹/² 𝑑² 𝑑𝜇² 𝑃�,� ≈ ‒ �𝑞² sin�𝑞𝑡𝐶��exp(‒𝜇𝑞)𝐽₀�𝑞𝑟�𝑑𝑞,𝛼𝐹₀𝑇(2 ‒ 𝜉²)𝐶� 4𝜋𝜌𝑄(𝜉)𝐶�² ∞ ₀ 𝑓(𝑥, 𝜇) = , 𝑔 = �𝑢�1‒𝑥₁² � + 𝑥₁�𝑥₁² ‒3��, 𝑢 = �1+ 𝑥₁² �¹/².3(𝑢 + 𝑥₁)¹/² 4𝑢⁵ 𝑔(𝑥₁) 𝜇⁵/² Vadim V. Surkov et al. 6 4. 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: (28) 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. 𝐸� = ‒ 𝐸� , 𝐸� = ,𝜑� 𝜇 𝑑𝑔 𝑑𝑥₁ = �𝑢𝑥₁�𝑥₁² ‒ 3� ‒ 𝑥₁⁴ + 6𝑥₁² ‒ 1�.𝑑𝑔 𝑑𝑥₁ 15(𝑢 + 𝑥₁)¹/² 8𝑢⁷ 𝐸� ~ .𝛼𝐶��𝐹₀𝑇(2 ‒ 𝜉²)𝐶� 8𝑟¹/² 𝑟₀³/² 𝑄(𝜉)𝜌𝐶�² 𝜑 = 𝜑�𝑔(𝑥1), 𝜑� = . 𝛼𝐶��𝐹₀𝑇(2 ‒ 𝜉²)𝐶� 8𝜋𝑟¹/² 𝜇⁵/² 𝜌𝐶�² 𝑄(𝜉) 7 Seismoelectric effect in Lamb’s problem Figure 1. A model calculation of (a) the dimensionless telluric potential and (b) the radial component of the electric field measured near ground surface as functions of the dimensionless coordinate. The seismoelectric effect is resulted from Rayleigh wave radiated by pulsed seismic source. The reference frame is used which moves at Rayleigh wave velocity. a) b) 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]: (32) For the actual solids the value of 𝜉 varies in the range from 0.874 to 0.955. Considering the upper layer of sedimentary rocks, the typical values of the Poisson coefficient are the following: 0.27 for large fragmental ground, 0.30 for sand and sandy loam, 0.35 for loamy soil and 0.42 for clay [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) we obtain that 𝜉 = 0.928, 𝛾 = 1.870 and 𝑄 = ‒1.314. The electrokinetic coefficient 𝐶�� can vary within a very wide range depending on the porosity, permeability, moisture saturation, concentration of mineral salts in the pore fluid and 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. 5. 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 (35) 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 (𝑡 > 𝑇) (36) 𝐸� = ‒ 𝐸� �𝐺�𝜏� ,𝑑𝜏 𝑇 𝑑𝑔 𝑑𝑥₁ ∞ ₀ 𝑓₁(𝑥, 𝜇) = ��(𝜇² + 𝑥²)¹/² ‒ 𝑥�¹/²‒�(𝜇² + 𝑥�² )¹/² ‒𝑥��¹/²�, 𝜕² 𝜕𝜇² 2 𝐶� 𝛾 = � �¹/².2(1 ‒ 𝜈) 1 ‒ 2𝜈 ∞ ₀𝑓₁(𝑥, 𝜇) =� 𝐺�𝜏�𝑑𝜏, 𝜕² 𝜕𝜇² �(𝜇² + 𝑥�² )¹/² ‒ 𝑥��¹/² (𝜇² + 𝑥�² )¹/² 𝜉 = .0.87 + 1.12 𝜈 1 + 𝜈 Vadim V. Surkov et al. 8 where 𝑥� = 𝑟 ‒ 𝐶� (𝑡‒𝑇) = 𝑥 + 𝐶�𝑇. Taking derivatives of this expression with respect to 𝜇 and rearranging one can reduce equation (36) to the view: 𝑓₁(𝑥, 𝜇) =[𝑔₁(𝑥₁) ‒𝑔₁ (𝑥�/𝜇)]/(2𝐶�𝜇³/²), (37) where (38) By making calculations analogous to those given above, one can find the potential and the radial component of the electric field (39) (40) Here as before the peak values of the telluric potential and horizontal component of the electric field are of the order of 𝜑� and 𝐸� = 𝜑� / 𝜇, respectively. The dimensionless potential 𝜑/𝜑� and electric field 𝐸� /𝐸� are shown in Figures 2a and 2b as functions of the dimensionless distance 𝑥₁. The signals shown with lines 1 correspond to the value 𝑥� / 𝜇 = 0.1 and are similar to the ones displayed in Figures 1a and 1b originated from the pulse source. The qualitative difference between the electric signals produced by the different sources becomes evident when considering lines 2 and 3, which correspond to the values 𝑥� / 𝜇 = 3 and 10, respectively. The graphs in Figure 2 are more complicated in shape than the ones shown in Figure 1. The individual peaks of the electric signals shown with lines 3 form two groups. One group is concentrated near the anterior front of the seismic wave whereas the other one near the posterior front. In the time interval between these groups, the electric signals are small because of small value of the pore pressure gradient. 𝑔₁(𝑥₁) = .𝑥₁² ‒ 𝑢𝑥₁ ‒ 1 𝑢³(𝑢 ‒ 𝑥₁)¹/² 𝜑 = 𝜑�[𝑔₁(𝑥₁) ‒𝑔₁(𝑥�/𝜇)], 𝜑� = ,𝛼𝐶��𝐹₀(2 ‒ 𝜉²) 16𝜋𝑟¹/² 𝜇³/² 𝜌𝐶�² 𝑄(𝜉) 𝐸� = 𝐸�[𝑔₂(𝑥₁) ‒𝑔₂(𝑥�/𝜇)], 𝑔₂(𝑥₁) = .3[𝑥₁(𝑥₁ ‒ 𝑢) (2𝑥₁ + 𝑢)]+ 𝑢‒2𝑥₁] 2𝑢⁵(𝑢 ‒ 𝑥₁)¹/² 9 Seismoelectric effect in Lamb’s problem Figure 2. The same as in Figures 1a and 1b but for the case of the seismic source describing by a step function. Lines 1, 2 and 3 correspond to the values 𝑥� /𝜇 = 0.1, 3, and 10, respectively. a) b) 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 𝐺(𝑡) = exp(‒𝑡/𝑡�)sin(2𝜋𝑓𝑡), (41) where 𝑓 is the seismic wave frequency and 𝑡� denotes the characteristic decay time. Substituting equation (41) for 𝐺(𝑡) into equation (35) and performing the numerical integration, one can find the electric field. A model calculation of the normalized radial component of the electric field and the source function 𝐺(𝑡) are shown in Figure 3a,b with solid and dotted lines, respectively. In making these plots we have used the following numerical values of the parameters 𝑓 = 0.5, 0.8 Hz; 𝑡� = 10 s; 𝐶� = 3 km/s; 𝑟₀ = 5 km and 𝑧 = 0. The zero- time reference in Figure 3a,b corresponds to the time 𝑡 = (𝑟 ‒ 𝑟₀)/𝐶� of seismic wave arrival at the observation point. 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]. 6. 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 Figure 3. Numerical modeling of the co-seismic signal caused by quasi-oscillatory seismic wave; (a) 𝑓 = 0.5 Hz, (b) 𝑓 = 0.8 Hz. The normalized radial component of electric field and the source function 𝐺(𝑡) are shown with solid and dotted lines, respectively. a) b) 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. Acknowledgements. This study was partly supported by the grant No. 18-05-00108 from RFBR. We appreciate useful comments of all Reviewers. References Anisimov, S.V., M.B. Gokhberg, E.A. Ivanov, M.V. Pedanov, N.N. Rusakov, V.A. Troitskaya, and V.E. Goncharov (1985). Short period vibrations of the Earth’s electromagnetic field due to industrial explosion, Rep. USSR Acad. Sci. 281 (3), 556-559. Balasco M., V. Lapenna, G. Romano, A. Siniscalchi, T.A. Stabile, and L. Telesca (2014). Electric and magnetic field changes observed during a seismic swarm in Pollino Area (Southern Italy), Bull. Seism. Soc. America 104 (3), 1289-1298, doi:10.1785/0120130183. Brekhovskikh, L.M., (1980), Waves in Layered Media, Second Edition, Applied Mathematics and Mechanics 6, Translated by R. T. Beyer, Academic press, New York. Eleman, F. (1965). The response of magnetic instruments to earthquake waves, J. Geomagn. Geoelectr. 18, 43-72. Frenkel, Ya.I. (1944). On the theory of seismic and seismoelectric phenomena in a moist soil, J. Phys. (USSR) 8, 230-241. Gao, Y., D. Wang, C. Yao, W. Guan, H. Hu, W. Zhang, P. Tong, and Q. Yang (2019). Simulation of seismoelectric waves 11 Seismoelectric effect in Lamb’s problem using finite-difference frequency-domain method: 2D SHTE mode, Geophys. J. Int. 216 (1), 414-438. Gorbachev, L.P., and V.V. Surkov (1987). Perturbation of an external magnetic field by a Rayleigh surface wave, Magnetohydrodynamics 23, 117-125. Grandshteyn I.S., and I.M. Ryzhik (2007). Table of integrals, series, and products, Seventh edition, Ed. by A. Jeffrey and D. Zwillinger, Elsevier, Academic press, p. 743. Guglielmi, A.V., (1986). Excitation of electromagnetic field oscillations by elastic waves in conducting body, Geomagn. Aeron. 27, 467-470. Han, Q., and Z. Wang (2001). Time-domain simulation of SH-wave-induced electromagnetic field in heterogeneous porous media: A fast finite element algorithm, Geophysics 66 (2), 448-461. Honkura, Y., M. Matsushima, N. Oshiman, and et al. (2002). Small electric and magnetic signals observed before the arrival of seismic wave. Earth Planets Space 54 (12), e9-e12. Hu, H., and Y. Gao (2011). Electromagnetic field generated by a finite fault due to electrokinetic effect, J. Geophys. Res. 116, B08302, doi:10.1029/2010JB007958. Iyemori, T., T. Kamei, Y. Tanaka, M. Takeda, T. Hashimoto, T. Araki, T. Okamoto, K. Watanabe, N. Sumitomo and N. Oshiman (1996). Co-seismic geomagnetic variations observed at the 1995 Hyogoken–Nanbu earthquake, J. Geomagn. Geoelectr. 48, 1059-1070. Ivanov, A.G. (1940). Seismoelectric effect of the second kind. Proc. USSR Acad. Sci. Ser. Geogr. Geophys. 4, 699-726. Jouniaux, L., M.L. Bernard, M. Zamora, and J.P. Pozzi (2000). Streaming potential in volcanic rocks from Mount Pele´e, J. Geophys. Res. 105, 8391-8401. Knopoff, E.L., (1955). The interaction between elastic waves motion and a magnetic field in electrical conductors, J. Geophys. Res. 60, 617-629. Lamb, H. (1931), The dynamical theory of sound, 2nd edition, London, Edward Arnold & Co. Landau, L.D., and E.M. Lifshitz (1970). Theory of Elasticity, Volume 7 of A Course of Theoretical Physics, Pergamon Press, The second English edition, 1970, 165 p. Martner, S.T., and N.R. Sparks (1959). The electroseismic effect, Geophysics 24, 297-308. Mavko, G., T. Mukerji, and T. Dvorkin (2014). The rock physics handbook: tools for seismic analysis of porous media, 2nd edn. Cambridge University Press, Cambridge Mogi, T., Y. Tanaka, D.S. Widarto, E.M. Arsadi, N.T. Puspito, T. Nagao, W. Kanda, and S. Uyeda (2000). Geoelectric potential difference monitoring in southern Sumatra, Indonesia - Co-seismic change, Earth Planets Space 52, 245-252. Monachesi, L.B., F.I. Zyserman, and L. Jouniaux (2018). An analytical solution to assess the SH seismoelectric response of vadose zone. Geophys. J. Int. 213 (3), 1999-2019. Nagao, T., Y. Orihara, T. Yamaguchi, I. Takahashi, K. Hattori, Y. Noda, K. Sayanagi and S. Uyeda (2000). Co-seismic geoelectric potential changes observed in Japan, Geophys. Res. Lett. 27, 1535-1538. Nikitin I.S., A.V. Philimonov, and V.L. Yakushev (2013). Rayleigh wave propagation at the oblique impact of a meteorite on the Earth’s surface and their impact on buildings and structures, Computer research and modeling, 5, No. 6, 981-992 (in Russian). Petrashen, G.I., G.I. Marchuk, and K.I. Ogurtsov (1950). On Lamb problem for the case of half space, Memoir of Leningrad State University, Mathematical series, issue 21, No. 55, 71-118 (in Russian). Pride, S. (1994). Governing equations for the coupled electromagnetics and acoustics of porous media, Phys. Rev. B 50, 15678-15696. Pride S.R., and S. Garambois (2005). Electroseismic wave theory of Frenkel and more recent developments, J. Eng. Mech. 131, 898-907. Romano, G., M. Balasco, A. Siniscalchi, A.E. Pastoressa, V. Lapenna (2018). Robust analysis for the characterization of the seismo-electromagnetic signals observed in Southern Italy, Ann. Geophys. 61, Doi: 10.4401/ag-7811. Sorochan, E. A. (1985). Handbook of the designer. Bases, foundations and underground structures, 480 p. (in Russian). Surkov, V.V., (1997). The electromagnetic precursor of a seismic wave, Geomagn. Aeron. 37, 792-796. Surkov, V.V., and M. Hayakawa (2014). Ultra and Extremely Low Frequency Electromagnetic Fields, Springer Geophysics. Surkov, V.V., and V.A. Pilipenko (2015). Estimate of ULF electromagnetic noise caused by a fluid flow during seismic or volcano activity, Ann. Geophys. 58(6), S0655, https://doi.org/10.4401/ag-6767. Surkov, V.V., V.A. Pilipenko, and A.K. Sinha (2018). Possible mechanisms of co-seismic electromagnetic effect, Acta Vadim V. Surkov et al. 12 http://link.springer.com/bookseries/10173 http://link.springer.com/bookseries/10173 http://link.springer.com/bookseries/10173 Geodaetica et Geophysica 53, https://doi.org/10.1007/s40328-018-0211-6 2018. Teodorovich, G.I. (1958). The doctrine of the formation of sedimentary rocks, Leningrad, 572 p. (in Russian). Vinogradov, J., M.Z. Jaafar and M.D. Jackson (2010). Measurement of streaming potential coupling coefficient in sandstones saturated with natural and artificial brines at high salinity, J. Geophys. Res., Solid Earth 115, No. B12, DOI:10.1029/2010JB007593. *CO R R E S PO N D I N G A U T H O R: Vadim Surkov, Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation, Russian Academy of Sciences, Moscow, Troitsk, Russia and Institute of Physics of the Earth, Russian Academy of Science, Moscow, Russia; e-mail: surkovvadim@yandex.ru © 2020 the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved 13 Seismoelectric effect in Lamb’s problem https://doi.org/10.1007/s40328-018-0211-6%202018 mailto:surkovvadim@yandex.ru