Semi-analytical solutions of seismo-electromagnetic signals arising from the motional induction in 3-D multi-layered media: part I—theoretical formulations

Taking into account the motional induction effect, in which the Earth crust that has finite electrical conductivity vibrates in the ambient geomagnetic field resulting in motionally induced electric current, we derive semi-analytical solutions of seismo-electromagnetic signals generated by an earthquake source in 3-D multi-layered media, which consists of an air half-space and multiple solid layers. First, both the elastic and electromagnetic (EM) wave-fields involved in the governing equations, which have the form of Maxwell’s equations coupled with elastodynamic equations, are expanded by a set of vector basis functions in cylindrical coordinate system. Then, we reorganize the transformed governing equations expressed by expansion coefficients and obtain corresponding first-order ordinary differential equations for the wave-fields in air and solid media. The expansion of the motionally induced electric current and the reorganization of Maxwell’s equations are the most important part, and also the most complicated and tedious part of this work. Thereafter, we solve the first-order ordinary differential equations through the Luco–Apsel–Chen generalized reflection and transmission method gaining solutions of the expansion coefficients. Finally, we obtain the frequency–space-domain semi-analytical solutions written as integrations of corresponding expansion coefficients over wavenumber domain, which can be numerically calculated by the discrete wavenumber method. The time-domain solutions can be achieved by further applying the discrete inverse Fourier transform. To have a numerical stability at any high frequency, we adopt the analytical regularization approach in the derivation process by introducing two artificial interfaces with infinitely small distance from the source. On the basis of the semi-analytical solutions, we can tell that only EM fields of TM mode (in which magnetic fields are transversely polarized) will be induced by SH waves, whereas EM fields of both TE mode (in which magnetic fields are transversely polarized) and TM mode will be induced by P and SV waves. The derived semi-analytical solutions can be used to calculate seismo-electromagnetic signals either below or above the free surface.


Introduction
Temporal variation of electromagnetic (EM) fields observed before and during earthquake events has been reported by many researchers (e.g., Belov et al. 1974;Gokhberg et al. 1989;Farser-Smith et al. 1990;Takeuchi Sun et al. Earth, Planets and Space (2021) 73:20 et al. 1997;Honkura et al. 2000Honkura et al. , 2004Nagao et al. 2000;Karakelian et al. 2002;Matsushima et al. 2002;Tang et al. 2010;Fujinawa et al. 2011;Han et al. 2011Han et al. , 2014Han et al. , 2015Han et al. , 2016Han et al. , 2017Huang 2011a, b;Hattori et al. 2013;Xu et al. 2013;Fujinawa and Noda 2015). These EM signals have attracted many attentions from geophysical scientists, because they are obviously connected to earthquake and provide a new perspective to detect and study development of earthquake and occurrence processes, as well as the precursory anomaly of earthquake. Hence, many theoretical and experimental efforts have been made to reveal the physical mechanism of the earthquakerelated EM signals. After decades of research, geophysical scientists have come up with several possible physical models that consider different seismo-electromagnetic coupling mechanisms to explain the earthquake-related EM signals. Nowadays, four kinds of mechanisms, i.e., electrokinetic effect (e.g., Mizutani et al. 1976;Pride 1994), piezoelectric effect (Ogawa and Utada 2000;Huang 2002), piezomagnetic effect (e.g., Okubo et al. 2011;Yamazaki 2016), and the motional induction effect (e.g., Gershenzon et al. 1993;Honkura et al. 2000Honkura et al. , 2004, are widely accepted and frequently used for analyzing the earthquake-related EM signals. These four mechanisms are based on relatively independent physical hypothesises, and in fact, it is hard to explain all observed EM signals by just one of the mechanisms. The realistic earthquake-related EM signals may be governed by part or all of the four mechanisms, and hence, it is worthwhile to conduct extensive investigations to figure out the seismoelectromagnetic conversion characteristics of the four different mechanisms. For example, a technique extending the Luco-Apsel-Chen (LAC) generalized reflection and transmission (GRT) method (Luco and Apsel 1983;Chen 1993), which is a well-known and very useful method for the study of synthetic seismogram (Martin and Thompson 1997;Chen 2007;Ge and Chen 2008), was introduced to simulate seismic waves and EM signals by considering the electrokinetic effect (Ren et al. 2007(Ren et al. , 2010a. Later, this technique was adopted to study the EM signals generated by earthquake sources (Ren 2009;Ren et al. 2012Ren et al. , 2015Ren et al. , 2016aRen et al. , b, 2020Zhang et al. 2013;Huang et al. 2015;Sun et al. 2019).
Motional induction effect is also known as seismodynamo effect (Honkura et al. 2000(Honkura et al. , 2009 or crustal dynamo effect (Iyemori et al. 1996). It is based on the physical idea that the seismic ground motion of Earth's conducting crust immersed in the ambient geomagnetic field can yield an electromotive force and a motional induction current, which causes temporal variations in EM fields. Gershenzon et al. (1993) investigated several seismo-electromagnetic coupling effects using a dipole model for the source of the earthquake-related EM signals to estimate the EM signals caused by motional induction, piezomagmetic and electrokinetic effects. Iyemori et al. (1996) adopted the motional induction effect to discuss coseismic EM signals of the 1995 Hyogoken-Nanbu M7.2 earthquake by considering an electrically conducting slab moving in geomagnetic field. Matsushima et al. (2002) recorded the coseismic EM signals of the 1999 M7.4 İzmit earthquake and its aftershock, and claimed that the seismo-dynamo effect is a plausible mechanism for the recorded EM signals. Ujihara et al. (2004) conducted magnetotelluric observations to record EM data for the aftershock of the 26 May 2003 M7.1 off Miyagi Perfecture earthquake, and they concluded that the seismo-dynamo effect is responsible for the EM signals associated with seismic waves. Tang et al. (2010) also proposed that the coseismic EM signals of the aftershocks of 2008 Wenchuan Ms8.0 earthquake were possibly caused by motional induction effect.
To get a fundamental insight into the motional induction effect, quantitative investigations were also conducted by researchers. Yamazaki (2012) developed a solution for motional induction effect in 3-D horizontally layered model with seismic waves approximated as plane waves. He showed that the amplitudes of magnetic field arising from a Rayleigh wave with displacement of 10 cm and a period of 30 s were as large as 0.1 nT. Gao et al. (2014) provided analytical solutions for motional induction effect in a homogeneous full-space model. Their analytical solutions also showed that the coseismic electric and magnetic signals corresponding to the seismic waves with the acceleration of 0.1 m/s 2 can be as large as 1 mV/km and 0.1 nT, respectively. Later, Gao et al. (2019) studied the EM responses to an earthquake source due to the motional induction effect in a 2-D horizontally layered model. These quantitative studies suggested that the motional induction effect possibly plays an important role in generated EM signals during earthquake. However, these studies have some limitations. The solutions of Yamazaki (2012) cannot apply for the near-field case, since the approximation of seismic plane wave was used. The analytical solutions (Gao et al. 2014) are only for a homogeneous full-space. The 2-D solutions of Gao et al. (2019) only account for a special situation, in which the fault plane is perpendicular to the plane of seismic wave propagation. Besides, Gao et al. (2019) claimed that it was difficult to derive the 3-D solutions by methods like the LAC GRT method, and that double space-to-wavenumber Fourier transforms must be applied to derive the 3-D solutions. However, they were wrong about this point.
In this study, we consider a 3-D multi-layered model, which is widely used in the study of seismology, and use the LAC GRT method to derive the semi-analytical solutions of seismo-electromagnetic signals generated by an earthquake source due to the motional induction effect. The mathematical derivations are shown step by step to assure the credibility. The remainder of this paper is organized as follows. First, we summarize the governing equations that have the form of Maxwell's equations coupled with elastodynamic equations, as well as boundary conditions for seismic waves and EM signals in a multilayered model. Second, we give a brief introduction to a set of vector basis functions, which is expressed in cylindrical coordinates, and use it to obtain transformed governing equations. In this process, deriving the expansion coefficients of the cross product of seismic displacement to ambient geomagnetic field is complicated but a key step. Later, we show the mathematical derivation process for obtaining the general solutions of the transformed governing equations in the frequency domain. Finally, we utilize the LAC GRT method to determine the source terms and wave-amplitude vectors, and insert them into the general solutions obtaining the semi-analytical solutions in frequency-space domain. The time-spacedomain solutions can be achieved by further applying the discrete inverse Fourier transform.

Governing equations and boundary conditions for multi-layered media
The multi-layered media concerned in the study are illustrated in Fig. 1. There are N homogeneous solid layers below an air half-space, and each solid layer is bounded by horizontally flat interfaces, z = z (j) and z = z (j+1) . The interface with the depth of z = z (1) is a free surface. The bottom layer (N-th solid layer) is a half-space and thus has a lower interface z = z (N +1) = +∞ . A double couple point source is embedded in the sth solid layer and located at (0, 0, z s ) . To eliminate the exponential growth factors in the source term (see Source term section), the analytical regularization approach (Chen 1999;Ren et al. 2010b) is applied by inserting two artificial interfaces z = z (s−) = z s − △ z and z = z (s+) = z s + △ z , where  z is positive and infinitely small quantity. When △ z approaches to zero (i.e., △ z → 0 ), the sth solid layer actually will be divided into two layers, z (s) ≤ z ≤ z (s−) and z (s+) ≤ z ≤ z (s+1) .
The elastic and electrical properties of each solid layer are described by Lamé constant (j) , shear modulus G (j) , density ρ (j) , magnetic permeability µ (j) , electrical permittivity ε (j) , and conductivity σ (j) . The two parameters µ (j) and ε (j) can be related to µ 0 (the vacuum magnetic permeability) and ε 0 (the vacuum electrical permittivity) by µ (j) = µ (j) r and ε (j) r represent the relative magnetic permeability and relative electrical permittivity, respectively. We will also derive the semi-analytical solutions of EM fields in the near-free-surface air (e.g., with altitude less than several kilometers), which is far away from ionosphere, and therefore presumably can be treated as non-conductive, i.e., σ (0) = 0 . For the air medium, we assume µ (0) = µ 0 and ε (0) = ε 0 .
Assuming a time dependence of e −iωt , one can find the governing equations of seismo-electromagnetic signals in the jth solid layer can be written as the following frequency-domain expressions: where j = 1, 2, . . . , N ; ε (j) = ε (j) + i ω σ (j) ; ω represents the circular frequency; u is the seismic displacement (meaning −iωu represents the seismic vibration velocity); F indicates the body force density; δ j,s is Kronecker delta function; τ is the traction acting on a horizontal plane below the air half-space; e z represents the unit vector along z-direction; E and H are the electric and magnetic fields, respectively; B a indicates the ambient geomagnetic field. The seismic waves are independent of EM signals. The second term in the right side of Eq.
(3) represents the induction electric current generated by the vibration of conductive media in geomagnetic field. It acts as the source of EM fields. This study considers time and spatial scales up to several minutes and several hundred kilometers, respectively. Although the ambient geomagnetic field B a varies over these scales, the amplitudes of these temporal-spatial variations are significantly smaller than the temporal-spatial average. (1) Therefore, temporal-spatial variations in B a should cause only minor changes on the resultant EM variations. For this reason, B a is considered as a uniform vector over the area and period of interest (Yamazaki 2012).
Equation (2) is derived from τ (j) = Ŵ (j) · e z , where Ŵ (j) is the stress tensor given by: where I is the identity matrix and the superscript 'T' represents the transpose of the matrix. The boundary conditions that should be taken into consideration in deriving the semi-analytical solutions are the reason of using τ instead of Ŵ (j) .
The detailed boundary conditions are as follows i. Traction free on the free surface, that is: ii. Continuities of seismic displacement and traction fields at subsurface interfaces, that is: iiii. Continuities of horizontal components of EM fields at each interface, that is: iv. Radiation condition at infinity ( z = ±∞ ), that is: where E (0) and H (0) indicate the EM fields in the air.

Transformed governing equations
In this section, we adopt a classic method that has been widely used as a pretreatment before deriving semi-analytical solutions of wave-fields in 3D multi-layered media (e.g., Aki and Richards 1980;Chen 1999;Ren et al. 2007Ren et al. , 2009Ren et al. , 2010aRen et al. , b, 2012. This pretreatment leads to transformed governing equations. The details are introduced below. , (j = 1, 2, . . . , N ); Sun et al. Earth, Planets and Space (2021) 73:20 The involved wave-field vectors are expanded using a set of vector basis functions: where Y m k (r, θ) = J m (kr)e imθ (m = 0, ±1, ±2, . . . ) is the cylindrical harmonic function. The orthogonality and completeness of this set of vector basis functions, which have been proven in previous studies (Aki and Richards 1980;Chen 1999), can be found in Appendix A.
An arbitrary vector, for example, C(r, θ, z) , can be expressed in terms of this set of function basis as: where C T ,m (z, k) , C S,m (z, k) , and C R,m (z, k) are the expansion coefficients given by: where the symbol * indicates the complex conjugate. (11)

Transformed governing equations of elastic waves
For the double couple point source concerned in this study, the body force density F can be mathematically expressed as: where r is the horizontal space variable, i.e., r = (r, θ) = (x, y) ; r s denotes the horizontal space variable of source point; z s indicates the source depth; δ is the delta function; M(ω) represents the spectrum of the source moment tensor. There is a trick in evaluating expansion coefficients of the source according to Eqs.
(13)-(15): if an arbitrary cylindrical coordinate system is chosen, the order 'm' will theoretically go to infinity. However, if a 'source-center' cylindrical coordinate system, in which the point source is located in z-axis, is chosen, the order 'm' will be limited as |m| ≤ 2 . Therefore, we choose a source-center coordinate system for 3D problem (see Fig. 1). Then, the expansion coefficients of the double couple point source can be written as: where: M xx (ω) , M xy (ω) , M xz (ω) , M yy (ω) , M yz (ω) , and M zz (ω) are the corresponding elements of the source moment tensor M(ω) . In this work, we take into account double couple point source representing a fault without tension (i.e., the displacement discontinuity across the fault is parallel to fault plane). For such kind of source, Eqs.
(3.21) and (3.23) in Chapter 3 of Aki and Richards (1980) tell us that the source moment tensor is a symmetric tensor (Chen 1999).
The elastic wave-fields involved in Eqs. (1) and (2) can be rewritten in terms of expansion coefficients. Using the orthogonality of the vector basis functions (see Appendix A) and noticing Eqs. (17) and (18), we find that the expansion coefficients of u and τ satisfy the following linear superposition: where: where: . We can see, the two expansion coefficients associated with SH waves, u (j) T ,m (z, k) and τ (j) T ,m (z, k) , are independent of those associated with P and SV waves, u R,m (z, k) . This fact leads to two uncoupled modes of elastic wave-fields; that is, SH mode governed by Eqs. (21)-(23) and PSV mode governed by Eqs. (24)-(28).

Expansion coefficients of the cross product u × B a
The induction electric current −iωσ u × B a acts as a source of EM fields. Determining the expansion coefficients of u × B a is the key of transforming Eq. (3) and solving EM fields. For 3-D case, the intersection angle between propagation direction of seismic wave and orientation of the ambient geomagnetic field varies for different azimuth angle. This makes the expansion of the vector u × B a and the derivation of EM field solutions in 3-D case more complicated and tedious than the 2-D case.
Now, let us define: Note that seismic displacement can be expressed by expansion coefficients as: and Eq. (11) can be rewritten as: where e r and e θ are the two horizontal unit vectors of cylindrical coordinate system.
Assuming the x-, y-and z-components of B a are B a x , B a y , and B a z , respectively, we can obtain the expression of B a in cylindrical coordinate system as: where, ξ = T , S, R ; and For the above nine integration formulas, we need to find out their explicit expressions one by one.
For Eq. (34), we can replace T m k (r, θ) , T m ′ k ′ (r, θ) , and B a with their expressions in cylindrical coordinate system. Utilizing Eqs. (31) and (32), we derive: Then, we substitute Eq. (43) into Eq. (34) and use the relation 1 2π 2π 0 e i(m−m ′ )θ dθ = δ m,m ′ to gain: Utilizing the following property of Bessel function we find: Combination of Eqs. (31) and (32) also gives: Substituting the above formula into Eq. (35) results in: Notice another property of Bessel function we gain: From Eqs. (31) and (32), we can also obtain: (52) (56) In a similar way, utilizing Eqs. (31)-(33), (37)-(42) and some properties of Bessel function, such as Eqs. (45), (49), (52), and (54), we gain: and Based on Eqs. (12) and (56)- (58), we find C(r, θ, z) can be divided into three parts, which are related with seismic displacement expansion coefficients of m ′ , m ′ + 1 , and m ′ − 1 orders, respectively. For these three parts, we use m to replace m ′ , m ′ + 1 and m ′ − 1 , respectively, to make the seismic displacement expansion coefficients have a uniform order; that is, u ξ ,m (z, k) ( ξ = T , S, R ). Since the expansion coefficients of u T ,m (z, k) and τ T ,m (z, k) (which correspond to SH waves) are not coupled with those of u S,m (z, k) , u R,m (z, k) , τ S,m (z, k) , and τ R,m (z, k) (which correspond to P and SV waves), the SH and PSV modes will be solved separately. Correspondingly, we separate the induction electric current induced by SH waves from those induced by P and SV waves. Finally, C(r, θ, z) can be divided into six parts as: Each part can be expressed in terms of expansion coefficients as: where, p = 0, +1, −1 ; ζ = SH , PSV ; and

Transformed governing equations of EM fields
Each part of the induction electric current, i.e., −iωσ C p,ζ (r, θ , z) , correspondingly will induce partial EM fields that contribute to the total EM fields. From this viewpoint, the total EM fields can also be divided into six parts and written as: Each part of EM fields can be expressed by expansion coefficients as: where, j = 0, 1, 2, . . . , N ; p = 0, +1, −1 ; ζ = SH , PSV . The EM fields in solid media (i.e., j = 1, 2, . . . , N ) satisfy: Rewriting the above two equations in terms of expansion coefficients and utilizing the orthogonality of the vector basis functions (see Appendix A), we obtain the general form of transformed governing equations of EM fields in solid (see Appendix B). Equations (170)-(172) involve the EM fields of TM mode (in which magnetic fields are transversely polarized), whereas Eqs. (173)-(175) involve the EM fields of TE mode (in which electric fields are transversely polarized).
Taking into account the specific expansion coefficients of C p,ζ (r, θ , z) , i.e., Eqs. (60)-(64), we find that SH waves will not generate EM fields of TE mode, that is: However, SH waves will generate EM fields of TM mode. The expansion coefficients of these TM-mode EM fields are actually determined by u T ,m (z, k) . Further considering the linear superposition of u T ,m (z, k) , i.e., Eq. (21), we obtain: SH waves with the case of p = +1, −1 determine another partial TM-mode EM fields, i.e., However, another two parts of induction electric current, −iωσ C ±1,PSV , will induce TM-mode EM fields. The expansion coefficients of these EM fields are related to u S,m (z, k) and u R,m (z, k) . Further considering the linear superposition of these two seismic-displacement expansion coefficients, i.e., Eq. (24), we obtain: We learn from Eqs. (170)-(172) that partial induction electric current −iωσ C 0,PSV does not generate EM fields T ,n (z, k).
Equations (173)-(175) tell us that PSV-mode seismic waves will also induce TE-mode EM fields because of the three non-zero expansion coefficients, C 0,PSV T ,m (z, k) , S,n (z, k).
PSV-mode seismic waves with the case of p = 0 determine partial TE-mode EM fields, i.e., fields' horizontal components at free surface. Using a similar way of deriving Eqs. (71)-(81), we can obtain the transformed governing equations of EM fields in air (see Appendix C).

Transformed governing equations written in matrix
SH waves are independent of P and SV waves for layered media considered in this study. Consequently, the transformed governing equations of displacement-stress-EM-wave-fields in solid media can be written as two sets of first-order ordinary differential equations associated with SH-mode and PSV-mode seismic waves, respectively. Although these two equation sets are different in dimension, they are identical in form: Since the air (above the free surface) is treated as nonconductive, i.e., σ (0) = 0 , we do not consider any induction electric current in the air. The generation of EM fields in air actually results from the continuities of EM where ζ = SH , PSV ; j = 1, 2, . . . , N.
Matrices ζ and ζ (z) are related with the eigendecomposition of matrix A ζ . The eigenvalues of A ζ actually have physical meanings: they are merely i times the vertical wavenumbers for down-going and up-going wave-fields (Aki and Richards 1980). In the following context, we use subscripts 'd' and 'u' to indicate the matrices (or vectors) associated with down-going and upgoing wave-fields, respectively.

Source term
After inserting the two artificial interfaces z = z (s−) = z s − △ z and z = z (s+) = z s + △ z , the source term s ζ is given by (Chen 1999;Ren et al. 2010b If the two artificial interfaces are excluded, the integration limits in Eq. (132), i.e., z (s−) and z (s+) , will be replaced by z (s) and z (s+1) , respectively. Besides, the two . That will make the source term contain exponential growth factors, such as e γ (s) s (z s −z (s) ) and e γ (s) s (z (s+1) −z s ) , leading to high-frequency instability problem (i.e., numerical overflow for high frequency) in the numerical calculation of the source term (Chen 1999;Ren et al. 2010b). The introduction of the two artificial interfaces z = z (s−) → z s − 0 and z = z (s+) → z s + 0 (as △ z → 0 ) eliminates the exponential growth factors and therefore helps to solve the highfrequency instability problem without sacrificing any computational efficiency (Chen 1999;Ren et al. 2010b  � ,
transmission due to the existence of interfaces (see Fig. 2). For the case of ζ = SH and PSV , they are 3 × 3 and 5 × 5 matrices, respectively, except that T SH (1) u and T PSV (1) u are 2 × 3 and 3 × 5 matrices, respectively.
The continuity boundary conditions, i.e., Eqs. (7) and (8), at subsurface interfaces help achieve the direct evaluation equations of the generalized reflection and transmission coefficients.

Fig. 2
Schematic diagram of the LAC GRT coefficients defined for different interfaces. a-c Coefficients defined for the real interfaces z = z (1) , z (2) , · · ·, z (N) , while (d) accounts for those defined for the two artificial interfaces z = z (s−) and z = z (s+) Since the Nth layer (bottom layer) is a half-space, the radiation boundary condition at z = z (N +1) = +∞ , i.e., Eq. (9), gives: On the free surface, the elastic wave-fields satisfy the traction-free condition, whereas the horizontal components of EM fields satisfy the continuity boundary condition. The combination of these two conditions leads to: For the case of ζ = SH , N SH 1 and N SH 2 , respectively, are 5 × 5 and 5 × 3 matrices given by: Now, the vector y ζ(j) (z) , which contains the expansion coefficients of elastic wave-fields and EM fields' horizontal components, can be computed for a receiver at any u , (j = 1, 2, . . . , s − 1).

The final solutions
After determining the source term and all the LAC GRT coefficients, we can calculate b 0, 1, . . . , s − 1 or j = s−, s+ or j = s + 1, s + 2, . . . , N ) through following equations:  (66) and (67), we obtain the semianalytical solutions of EM fields in the frequency-space domain. Similarly, the combination of Eqs. (21), (24), and (30) yields the semi-analytical solutions of seismic waves in the frequency-space domain. The wavenumber integration in Eqs. (30) and (67) can be numerically computed using the well-known discrete wavenumber method (Bouchon and Aki 1977;Bouchon 1981Bouchon , 2003. The final transformation back to the time domain can be performed by the discrete inverse Fourier transform.
It should be mentioned, if the source and receiver are located at close or same depths, it will be extremely time-consuming to numerically calculate the semi-analytical solutions without applying any other technique, because the convergence of wavenumber integration will become very slow. The peak-trough averaging method (Zhang et al. 2001(Zhang et al. , 2003, which is mathematically simple and easy to implement in practice, has been proven an effective and efficient method to overcome this slow convergence problem. Application of the peak-trough averaging method in our derived semi-analytical solutions will make it possible to efficiently calculate both seismic and EM signals generated by an earthquake with a shallow focus.

Conclusions
Using a set of vector basis functions in cylindrical coordinate system (Aki and Richards 1980;Chen 1993) to expand the involved wave-fields and solving the expansion coefficients through the LAC GRT method (Luco and Apsel 1983;Chen 1993Chen , 1999Ge and Chen 2008;Ren et al. 2007Ren et al. , 2010aRen et al. , b, 2012, we derive the semi-analytical solutions of seismo-electromagnetic signals arising from the motional induction in 3-D multi-layered media due to a double couple point source, as which an earthquake source can be treated under far-field condition (Aki and Richards 1980). The determination of the EM fields excited by the induction electric current in the 3-D case, which is the most complicated and tedious part of this work, is accomplished carefully and patiently.
The LAC GRT method adopted in the derivation has been proven numerically efficient and stable in previous studies (Martin and Thomson 1997;Chen 1999;Ge and Chen 2008). The analytical regularization approach (Chen 1999;Ren et al. 2010b), which introduces two artificial interfaces with distance △ z from the source and considers a limiting process of △ z → 0 , is adopted in an effort to solve the high-frequency instability problem that exists in the numerical computation. Therefore, the derived semi-analytical solutions should have advantage of higher efficiency and stability in numerical computation. Besides, we propose using the peak-trough averaging method (Zhang et al. 2001(Zhang et al. , 2003 to overcome another computational problem, that is, the slow convergence problem which occurs when the source and receiver are located at close or same depths. The derived semi-analytical solutions account for not only the seismo-electromagnetic signals in the solid but also those in the air. They indicate that SH waves only induce EM fields of TM mode, whereas P and SV waves induce EM fields of both TE and TM modes. On the basis of the semi-analytical solutions, expected characteristics of the seismo-electromagnetic signals arising from the motional induction in 3-D multi-layered media will be numerically investigated in a companion paper. and γ (0) em = k 2 − ω 2 µ 0 ε 0 with Re(γ (0) em ) > 0. The expansion coefficients of EM fields (in air) associated with P and SV waves are governed by: where n = 0, 1, 2 ; and where q = 0, 1 ; n = 0, 1, 2 ; and a TE(0) = 0 (γ (0) em ) 2 (iωµ 0 ) iωµ 0 0 .