Semi-analytical solutions of seismo-electromagnetic signals arising from the motional induction in 3-D multi-layered media: part II—numerical investigations

In this paper, numerical computations are carried out to investigate the seismo-electromagnetic signals arising from the motional induction effect due to an earthquake source embedded in 3-D multi-layered media. First, our numerical computation approach that combines discrete wavenumber method, peak-trough averaging method, and point source stacking method is introduced in detail. The peak-trough averaging method helps overcome the slow convergence problem, which occurs when the source–receiver depth difference is small, allowing us to consider any focus depth. The point source stacking method is used to deal with a finite fault. Later, an excellent agreement between our method and the curvilinear grid finite-difference method for the seismic wave solutions is found, which to a certain degree verifies the validity of our method. Thereafter, numerical computation results of an air–solid two-layer model show that both a receiver below and another one above the ground surface will record electromagnetic (EM) signals showing up at the same time as seismic waves, that is, the so-called coseismic EM signals. These results suggest that the in-air coseismic magnetic signals reported previously, which were recorded by induction coils hung on trees, can be explained by the motional induction effect or maybe other seismo-electromagnetic coupling mechanisms. Further investigations of wave-field snapshots and theoretical analysis suggest that the seismic-to-EM conversion caused by the motional induction effect will give birth to evanescent EM waves when seismic waves arrive at an interface with an incident angle greater than the critical angle θc = arcsin(Vsei/Vem), where Vsei and Vem are seismic wave velocity and EM wave velocity, respectively. The computed EM signals in air are found to have an excellent agreement with the theoretically predicted amplitude decay characteristic for a single frequency and single wavenumber. The evanescent EM waves originating from a subsurface interface of conductivity contrast will contribute to the coseismic EM signals. Thus, the conductivity at depth will affect the coseismic EM signals recorded nearby the ground surface. Finally, a fault rupture spreading to the ground surface, an unexamined case in previous numerical computations of seismo-electromagnetic signals, is considered. The computation results once again indicate the motional induction effect can contribute to the coseismic EM signals.


Introduction
In the companion paper (Sun et al. 2021), which is part I of current work, a set of systematic semi-analytical solutions of the seismo-electromagnetic signals arising from the motional induction in 3-D multi-layered media due to an earthquake source was obtained. In addition to the seismic and electromagnetic (EM) signals in the solid, EM signals in the air can also be calculated by the derived semi-analytical solutions. In the derivation process, Sun et al. (2021) utilized the linear superposition principle, applied the analytical regularization approach (Chen 1999;Ren et al. 2010a) by introducing two artificial interfaces infinitely close to the source, and adopted the Luco-Apsel-Chen (LAC) generalized reflection and transmission (GRT) method (Luco and Apsel 1983;Chen 1993Chen , 1999Ren et al. 2007Ren et al. , 2010bRen et al. , 2012Ren et al. , 2016aRen et al. , 2016bRen et al. , 2020. These measures bring higher efficiency and stability to the numerical computation of the semi-analytical solutions. Although the motional induction effect has been proposed as a possible generation mechanism of seismoelectromagnetic phenomena since nearly three decades ago (Gershenzon et al. 1993), quantitative investigations on this topic are very few. To our knowledge, only two works (Yamazaki 2012;Gao et al. 2019) provide solutions and conduct computations for seismo-electromagnetic signals arising from motional induction effect in multi-layered media. Considering a simple situation in which seismic waves are approximated by plane waves and the electrical conductivity of the Earth's crust has a stratified structure, Yamazaki (2012) derived analytical solutions of the EM responses and conducted numerical analysis. The numerical computations of Yamazaki (2012) aimed at some given periods (or frequencies) of seismic waves and did not require wavenumber integration. It should be noted that the solutions of Yamazaki (2012) are only applicable in the far-field condition because of the applied approximation of seismic plane waves. Considering 2-D case, Gao et al. (2019) used the global matrix method to derive the solutions of seismo-electromagnetic signals that are expressed as wavenumber integrations. They computed these wavenumber integrations to achieve the space-domain solutions. As shown by Eqs. (21) and (22) of Gao et al. (2019), 2-D solutions concerned in their work can be computed by inverse Fourier transform with respect to the wavenumber, because they applied a space-to-wavenumber Fourier transform to the involved wavefields. The 2-D solutions of Gao et al. (2019) are only applicable for a special case in which the fault plane is perpendicular to the wave propagation plane. For the 3-D case considered in this work, although the semi-analytical solutions (Sun et al. 2021) are also expressed as wavenumber integrations, they should be regarded as the summation of several inverse Hankel transforms (see Appendix A). Therefore, the 3-D semi-analytical solutions (Sun et al. 2021) cannot be computed by the inverse Fourier transform. Instead, they can be computed using the well-known discrete wavenumber method (Bouchon and Aki 1977;Bouchon 1981Bouchon , 2003. It should be mentioned, for both 2-D and 3-D cases, a slow convergence problem will occur in the numerical computation of the solutions when the source and receiver are located at close or same depths. This slow convergence problem has not been solved in the work of Gao et al. (2019).
The discrete wavenumber method, which introduces an infinite set of secondary sources of concentric rings distributed at equal radial intervals L p (spatial periodicity) to transform the wavenumber integration into summation (Bouchon and Aki 1977;Bouchon 1981Bouchon , 2003, has been applied in the numerical modelling studies of electrokinetic effect (e.g., Haartsen and Pride 1997;Garambois and Dietrich 2002;Ren et al. 2012). In these numerical computations of seismo-electromagnetic signals induced by electrokinetic effect, it was found using largest seismic wave velocity to determine the spatial periodicity L p is sufficient to guarantee the accuracy of numerical solutions. In this work, we also apply the discrete wavenumber method to numerically compute the semi-analytical solutions of Sun et al. (2021). We adopt the peak-trough averaging method (Zhang et al. 2001(Zhang et al. , 2003 to overcome the slow convergence problem. Besides, we utilize the point source stacking method (Olson and Apsel 1982;Ren et al. 2012) to deal with the case of a finite fault.
The remainder of this paper is organized as follows. First, we introduce our numerical computation approach enabling us to deal with any focus depth and to consider either a double couple point source or a finite fault. Second, we try to verify the proposed method by comparing with another numerical method for the seismic wave solutions. Last, we conduct numerical investigations that are different from those of Yamazaki (2012) and Gao et al. (2021). The coseismic characteristic of the EM signals, the generation of Keywords: Seismo-electromagnetic signals, Motional induction effect, Numerical computation, Coseismic EM signals, Evanescent EM waves evanescent EM waves on interfaces, and the effect of conductivity at depth on EM signals are presented and discussed.

Numerical computation approach
In the following context, we specify the detailed implementation of the discrete wavenumber method (Bouchon and Aki 1977;Bouchon 1981Bouchon , 2003, the peak-trough averaging method (Zhang et al. 2001(Zhang et al. , 2003, and the point source stacking method (Olson and Apsel 1982;Ren et al. 2012) in our numerical computation of seismoelectromagnetic signals arising from the motional induction due to a double couple point source or a finite fault embedded in a 3-D multi-layered media.

Discrete wavenumber method
According to Bouchon and Aki (1977) and Bouchon (1981Bouchon ( , 2003, the spatial periodicity of secondary sources L p determines the wavenumber sampling ∆k used in summation: Considering a cylindrical coordinates system with a source located at (0, 0, z s ) and a receiver located at (r, θ, z), the spatial periodicity L p should satisfy the condition that "pseudo" waves generated by the secondary sources do not enter the observational time window t ∈ [0, t max ], that is, there will be no interference from the nearest secondary source to the receiver's recording before the time of t = t max . This specific condition can be written as follows: where v max should be, in principle, the maximum wave velocity in the medium. Consequently, when considering the propagation of coupled seismic and EM waves, the value of v max theoretically should be the largest EM wave velocity, which implies a very large periodicity L p . Thus, according to Eq. (1), we should use a very small wavenumber sampling ∆k to ensure the correctness of the calculation results, which leads to unacceptably long computation time in practice. However, for the converted EM waves whose phase velocity is 1 √ µ(ε + iσ/ω) (where μ is magnetic permeability; ε is electrical permittivity; σ is conductivity; ω is the circular frequency), we have numerically verified that the computed synthetic signals are identical whether v max is based on the largest seismic or EM wave velocity.
(2) L p > 2r, The reason is that, although secondary sources also give birth to converted EM waves that propagate in both air and solid with EM wave velocity of the media, those converted EM waves contributed from secondary sources show a rapid amplitude decrease as the offsets between secondary sources and receiver increase. A spatial periodicity L p satisfying Eq. (2) with v max based on the largest seismic wave velocity can guarantee that the secondary source-receiver offsets are long enough for the attenuation of the converted EM waves from secondary sources. As a result, the amplitudes of EM waves generated by secondary sources are always negligible with respect to the EM waves generated by the principal source. Such a measure of using largest seismic wave velocity as v max to determine L p has been adopted in some previous studies that numerically calculated seismo-electromagnetic signals induced by electrokinetic effect (e.g., Haartsen and Pride 1997;Garambois and Dietrich 2002;Ren et al. 2012).
To perform the summation of secondary sources, one needs to remove the singularities of the integrands from the real k-axis. This can be achieved by giving the real circular frequency ω R a small imaginary part ω I . Hence, we use a complex circular frequency: when a time dependence of e −iωt is considered. In this way, the singularities are located in the first and third quadrants of the complex k plane on the straight line passing through the origin, eliminating the singularities on the real k-axis (Aki and Richards 1980). The value of ω I is usually chosen to be (Haartsen and Pride 1997): The introduction of complex frequencies leads to a smoothing effect on the spectrum, because it can relatively increase the strength of earlier-arriving signal. Besides, it can also reduce the interference generated from the nearest secondary source (Bouchon and Aki 1977). The influence of the imaginary part of the complex frequency can be removed from the final time-domain solutions using the inverse Fourier transformation with the kernel exp(−iω R t + ω I t).

Peak-trough averaging method
According to Sun et al. (2021), the semi-analytical solutions of seismic waves or EM signals due to an arbitrary point seismic source buried in a multi-layered media can be expressed as a summation of the products of the radiation patterns and the following type of oscillatory integrals: where r is the epicentral distance; z is the receiving point depth; z s indicates the source depth; k is the horizontal wavenumber; ω is the circular frequency; J n (kr) is the nth order Bessel function of the first kind; S(k, z, z s , ω ) represents the kernel function containing the attenuation factor. If the source-receiver depth difference becomes smaller, the attenuation factor contained in the kernel function will have weaker effect. Then, under the combined action of the kernel function and the Bessel function, the convergence of wavenumber integration becomes extremely slow. Thus, to obtain accurate computation results, one has to take a large upper limit of truncated integral, which leads to extremely low computation efficiency (Zhang et al. 2001(Zhang et al. , 2003. 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. When the wavenumber is greater than a certain critical value, the integral value oscillates around a certain value as the integral upper limit increases, and its envelope is a monotonous and smooth attenuation curve, which is characterized by slow convergence. The peak-trough averaging method can utilize a small number of peak and trough values to achieve the purpose of efficient and high-precision computation. The reason is that the peak-trough distribution of the slow convergence integral also appears as a monotonically decayed alternating sequence. Therefore, the peak and trough values of the integration curve can be extracted to form a slowconvergent sequence, whose convergence value (i.e., the integral value when the upper limit of the integral tends to infinity) can be obtained quickly and accurately using the repeated average method (Dahlquist and Björck 1974). The specific process of the peak-trough averaging method can be divided into the following 3 steps:

Determine the critical wave value k c
When k > k c , the envelope of the integral value appears as a monotonic attenuation curve, so we first determine the critical wave value k c , which can be given by the following empirical formula: (5) I n (r, z, z s , ω) = +∞ 0 S(k, z, z s , ω) · J n (kr)dk, where v min is the minimum velocity of the structure model; ω R is the real circular frequency; ω I is the imaginary circular frequency; α is an empirical coefficient, which is sufficiently large to ensure the decaying property of the integral at k > k c . As suggested by Zhang et al. (2001Zhang et al. ( , 2003, we use an empirical coefficient of 1.5 in our numerical computations.

Determine the peak and trough values
When k > k c , we first need to determine the initial peaks and troughs. Therefore, we record the positions of each adjacent three integral sampling points k i and corresponding function values (k i , G i ) (i = 1, 2, 3). If G 2 is greater or smaller than G 1 and G 3 , then (k 2 ,G 2 ) is an initial position of a peak or trough. If not, it proves that there is no peak or trough in this interval. In either case, move one point forward to form a new group of three successive integration sampling points and corresponding function values, and then repeat the above process, continuously increasing the wave number k by a fixed step ∆k, which is specified by Eq. (1), until enough peaks and troughs are obtained. The accuracy of the peak-trough averaging method depends on the accuracy of the function values of the peaks and troughs. However, due to the limitation of the step size, the accuracy of the initial peaks and troughs determined through the above screening process is often not accurate enough and must be further improved. Zhang et al. (2001Zhang et al. ( , 2003 proposed using the quadratic interpolation technique to find the accurate peaks and troughs within the allowable range of computation error. According to the above steps, enough initially determined peaks or troughs positions (k 2 , G 2 ) with their adjacent points (k 1 , G 1 ) and (k 3 , G 3 ) have been recorded. Although the location of the peak or trough is rough, it can be judged that the accurate peak or trough is located in the (k 1 , k 3 ) interval. Thus, the following quadratic interpolation polynomial can be constructed: Using the above method, a series of accurate peaks and troughs can be obtained.

Determine the integral convergence value
After finding a series of accurate peaks and troughs, and composing them into an oscillating slow convergence sequence, one can use the repeated average method (Dahlquist and Björck 1974) to obtain the convergence value, i.e., the required integral value.

The point source stacking method
It is known that, a fault can be considered as a double couple point source under the far-field condition. However, under near-field condition, a finite fault must be considered as a non-point source, because the pointsource approximation in this case will result in unacceptable errors of the wavefields computed (Ren et al. 2012). It is usually difficult to directly calculate the field generated by a finite fault. One solution to this problem is the point source stacking method (Olson and Apsel 1982;Ren et al. 2012). For a finite fault, the generated seismic or EM wavefields at one spatial point are the integral of Green's function over the area of the finite fault. In the numerical computation, integral is carried out by summation. Therefore, the finite fault can be discretized into a lot of cells. Each can be approximately treated as a point source located at the cell center when its area is sufficiently small. In this way, the finite fault can be represented by the serried distribution of numerous point sources in the fault area. Thus, the seismic or EM wavefields at one spatial point are the stacking of wavefields generated by each point source, which can be numerically computed through the semi-analytical solutions derived by Sun et al. (2021).
In theory, smaller cell size leads to more accurate numerical solutions. However, in actual numerical calculations, smaller cell size also means longer the calculation time. Therefore, we need to select a suitable cell size, which can not only reduce the computation time as much as possible but also satisfy the required precision. This aim can be achieved by the following steps (which take seismic wave solutions as an example): (i) Given a discretization scheme, for instance, ∆s = (L/M)(W/N) (where L and W indicate the length and width of the finite fault; M and N are positive integers; ∆s is the cell area), we compute the seismic wave solutions u (1) for the receiver point; (ii) Double the positive integers (2M → M and 2N → N) and recompute the seismic wave solutions u (2) ; (iii) Calculate the residual error e res through the following formula: (iv) Compare the residual error e res with the required precision p req . If e res < p req , the used cell size will be considered suitable and the computation results will be taken as the correct solutions. Otherwise, repeat the steps (ii)-(iv) until the correct solutions are obtained. The above steps should also be applied to the solutions of EM fields. A required precision of p req = 10 -6 is adopted in this study. In the following numerical computation, the North, East and downward directions are set to be x-, y-and z-directions, respectively.

Validation of the computed seismic waves
At present, there are still very few quantitative studies (Yamazaki 2012;Gao et al. 2019) of the seismo-electromagnetic signals arising from the motional induction effect. Since Yamazaki (2012) did not compute the timedomain full waveforms and Gao et al. (2019) only provided 2-D solutions, our work is the first one capable of computing full waveform of seismo-electromagnetic solutions in 3-D multi-layered media due to the motional induction. Therefore, we cannot find another method to compare both seismic waves and EM signals for this exact case. However, there are many existing methods, such as finite-difference method, that can compute seismic waves for the case concerned. The curvilinear grid finite-difference method (CGFDM) has been verified to be accurate for the modeling of seismic wave propagation in solid media (Zhang and Chen 2006;Zhang et al. 2012;Sun et al. 2016Sun et al. , 2018. Therefore, we compare our  numerical results with those obtained from the CGFDM by considering only seismic waves. A four-layer model (as illustrated in Fig. 1) is adopted. The four layers, from top to the bottom, consist of s 1 , s 2 , s 3 , and s 4 materials, respectively. Their elastic properties are listed in Table 1. A double couple point source is located at (0, 0, 10) km. The spectrum of the source moment tensor M(ω) can be expressed as where M xx , M xy , M xz , M yy , M yz , and M zz , are the moment tensor components and s(ω) is the source time function. The adopted source is assumed to represent a fault whose strike, dip, and rake angles are 30º, 72º, and 65º, respectively. The seismic moment is set to be 5.56 × 10 17 N·m, which corresponds to a M5.8 earthquake. Therefore, the moment tensor components are given by: M xx = − 2.67 × 10 17 N·m, M xy = 2.40 × 10 17 N·m, M xz = 1.41 × 10 17 N·m, M yy = -0.29 × 10 17 N·m, M yz = -3.89 × 10 17 N·m, and M zz = 2.96 × 10 17 N·m. In this study, the adopted source time function s(ω) is a Ricker wavelet which can be written as where ω is the circular frequency; f p and t 0 are the peak frequency and the time delay, respectively. For this fourlayer model, we set f p = 1 Hz and t 0 = 2 s.
Five receivers marked as triangles with specific numbers 1-5 (see Fig. 1) are considered. They are located on the ground surface z = 0 km and have the same y-coordinate y = 20 km. Their x-coordinates are x = 20, 30, 40, 50, 60 km. In Fig. 2, the seismic waves computed by our method (red lines) are compared with the CGFDM solutions (blue lines) for the five receivers with numbers 1-5. The three rows show the three components of seismic vibration velocity v x , v y , and v z . One can find a perfect agreement between the two sets of numerical results, which, to a certain extent, verify the correctness of the semi-analytical solutions derived by Sun et al. (2021) and the validity of the numerical computation approach proposed by us.
For current problem, both the seismic and EM solutions can be regarded as summation of several inverse Hankel transforms (see Appendix A). Due to this similarity of mathematical expressions, there is a high possibility that the discrete wavenumber method and peak-trough averaging method, which have been developed for calculating seismic waves, can also be applied to compute EM fields. A more complete and convincing validation that also takes into account the computed EM signals may be carried out in the future when another method capable of computing seismo-electromagnetic solutions for 3-D multi-layered media is available.

Numerical investigations
In the following section, we carry out numerical computations to investigate the characteristics of seismoelectromagnetic signals arising from motional induction. The air medium above the ground surface is assumed to be insulated and its magnetic permeability and electrical permittivity are set to be those values of the vacuum.
Attenuations which are often represented by quality factors exhibit for both P and S waves in the real Earth (Press 1964). Therefore, quality factors of P and S waves (i.e., Q p and Q s ) are considered in the following numerical  Fig. 2 Comparison of the synthetic v x , v y , and v z components for the five receivers marked as triangles with specific numbers 1-5 (see Fig. 1). The read lines indicate the results computed by our method, while the blue lines indicate the numerical solutions obtained by the CGFDM computation. Correspondingly, the slownesses of P and S waves (s p and s s ) are calculated by following formulas: where ρ, λ and G represent density, Lamé constant and shear modulus, respectively.

Synthetic recordings of an air-solid two-layer model
To facilitate analysis, we first consider a simple model, an air-solid two-layer model, which consists of an air upper half-space and a solid lower half-space. The solid lower half-space is made up of s 5 material, whose properties are listed in Table 2. The ground surface is chosen to be the plane z = 0 km. The source is a double couple point source located at (0, 0, 40) km. It is assumed to represent a fault with strike, dip, and rake angles of 0°, 72°, and 65°. Its seismic moment is 2.76 × 10 17 N·m, which corresponds to a M5.6 earthquake. The ambient geomagnetic field has an intensity of 5 × 10 -5 T. Its inclination and declination angles are 45° and 10°, respectively. The source time function is a Ricker wavelet with a peak frequency of 0.8 Hz and a time delay of 2.5 s. Seismo-electromagnetic signals are checked for two receivers nearby the ground surface. One is buried 0.1 m underground, while the other one is 2 m above the ground surface. Both of the two receivers have the same horizontal coordinates of x = 95.76 km and y = 80.35 km. Figure 3 shows the seismic vibration velocity (v x , v y and v z ), electric field (E x , E y and E z ), and magnetic induction intensity (B x , B y and B z ) computed for the (11) underground receiver (yellow lines), as well as EM signals (E x , E y , E z , B x , B y , and B z ) calculated for the other receiver in the air (thin blue lines). In exception to the E z components, the other five EM components (E x , E y , B x , B y , and B z ) are continuous across an interface. Besides, the two receivers are close enough, since their distance is negligible compared with the seismic wavelength. As a result, an overlap of the yellow and blue lines occurs for those five EM components (E x , E y , B x , B y , and B z ). Seismic arrivals of P and S waves, which start to show up at ~ 24.5 s and ~ 42 s, respectively, are obvious in both seismic and EM signals. Therefore, the EM signals displayed in Fig. 3 show up simultaneously with seismic arrivals. They are the so-called coseismic EM signals. This numerical result confirms that the motional induction effect is one of the possible generation mechanisms of coseismic EM signals associated with natural earthquakes (Iyemori et al. 1996;Honkura et al. 2000Honkura et al. , 2004Matsushima et al. 2002;Ujihara et al. 2004;Tang et al. 2010). Since the blue lines (in Fig. 3) correspond to the receiver above the ground surface, our numerical results suggest that coseismic EM signals can be recorded by receivers in the air. Actually, such a phenomenon was already reported in field observation by Ujihara et al. (2004). They set up induction coils hung in the air between trees to record magnetic signals. The observed magnetic field after band-pass filtering operation also exhibited the coseismic characteristic. Therefore, our computational result indicates that motional induction effect could be one of the possible generation mechanisms for the coseismic magnetic signals recorded in the Earth's near-surface air ). Other possible generation mechanisms include the electrokinetic effect, the piezomagnetic effect and so on. Further efforts to validate this idea could be made in the future.

Wave-field snapshots
Still using the air-solid two-layer model, we investigate the propagation of the seismic and EM wavefields by computing snapshots at different times for a rectangular area that has an epicentral (radial) distance range from 0 to 100 km, an azimuthal angle of 40°, and a depth range from − 5 to 15 km. The wave-field snapshots are determined at 401 × 81 receiver positions. The receiver spacing is 0.25 km in either epicentral (radial) direction or vertical direction. Figure 4 displays the snapshots of seismic vibration velocity v x component at times of 10.23, 14.23, 18.23, Fig. 4 Snapshots of seismic vibration velocity v x component at times of 10.23, 14.23, 18.23, 22.23, and 26.23 s. They are displayed for a rectangular area with an epicentral distance range from 0 to 100 km, an azimuthal angle of 40°, and a depth range from -5 to 15 km 22.23, and 26.23 s. In the top snapshot (i.e., at the time of t = 10.23 s), the direct S wave is still propagating upward, whereas the direct P wave has arrived at the ground surface (z = 0 km) giving birth to the reflected PP and PSV waves that start to propagate downward. At the time of t = 14.23 s (the second snapshot), the direct S wave just reaches the ground surface and the reflected PP and PSV waves are separated from each other, since they have propagated downward for a while. The third snapshot (t = 18.23 s) exhibits the SS and SVP waves reflected downward off the ground surface in addition to the direct P and S waves and reflected PP and PSV waves. The signal strength of the reflected PP wave becomes much weaker in comparison with other waves. In the bottom two snapshots computed at times of t = 22.23 s and t = 26.23 s, all the waves further propagate outward as time elapses.
Under the ambient geomagnetic field B a , seismic vibration in a conductive medium gives birth to the induction electric current σ v × B a , where σ and v are the material conductivity and the seismic vibration velocity, respectively. This induction electric current causes accumulation and depletion of electric charges that naturally result in an electric field. Under this electric field, electric charges has to move forming the conductive current σ E . In a homogeneous medium, the induction electric current and the conductive current will be counterbalanced by each other and an equilibrium will be reached eventually. The curl electric field will further generate magnetic field. Therefore, EM fields accompanying seismic waves will be induced in a conductive medium. They are local responses to the seismic waves passing by and hence can be called localized EM fields. Figures 5 and 6 show the snapshots of electric field E x component and magnetic field B x component, respectively, for the same times and the same area as Fiugre 4. For the area below the ground surface, one can easily identified part of the localized electric field, which accompanies the S waves (including the direct S wave, the reflected PSV wave and the reflected SS wave), as well as the localized magnetic field, which accompanies all seismic waves. The localized electric field accompanying the P waves (including the direct P wave, the reflected PP wave and the reflected SVP wave) actually also exists but its signal strength is much weaker than those associated with the S waves. The cause of this result is that of the solid half-space (s 5 material) has a conductivity of 0.001 S/m. For such a conductivity value, Gao et al. (2014) showed that the localized electric field of P waves could be more than two orders of magnitude weaker than that of S waves at the frequency of 1 Hz. Hence, the signal strength characteristic of the localized electric field in Fig. 5 is consistent with the analysis of Gao et al. (2014). Besides, Gao et al. (2014) also showed that the P and S waves' capacity of inducing localized electric field generally become closer when the conductivity increases from 0.0001 to 1 S/m for a frequency of 0.1, 1, 10, or 100 Hz. We also computed the snapshots by changing the conductivity of the solid half-space to 0.1 S/m. For that case, the localized electric field associated with both P and S waves can be clearly identified. The snapshots of that case are not displayed here to save space.
The localized EM fields exist in the solid but is absent in the air. Thus, they cannot be continuous across the ground surface. To satisfy the EM boundary condition, i.e., the continuity requirement of the total EM tangential components at an interface, corresponding non-localized EM fields are then generated. As a result, EM fields can also be observed in the air. Those in-air EM fields decay rapidly when moving away from the ground surface (see Figs. 5 and 6).
If we use k to indicate the horizontal wavenumber of seismic waves arriving at the ground surface, the localized EM fields will have the same horizontal wavenumber, which is given by k = sin θ(ω V sei ) , where V sei is the seismic wave velocity, and θ is the incident angle of the seismic waves arriving at the ground surface. Thereafter, as a result of the continuity boundary condition of EM fields at the ground surface, the horizontal wavenumber of non-localized EM fields in the air should also be equal to k. This means the horizontal wavenumber also satisfies the equation k 2 − (γ em ) 2 = (ω/V em ) 2 , where V em is the speed of light when we consider the EM fields in the air and γ em = k 2 − (ω/V em ) 2 is i times the vertical wavenumber of the non-localized EM fields. There is no down-going EM fields in the top layer (i.e., in the air); thus, a depth-dependent factor exp(γ em z rcv ) (where z rcv indicates the depth of an in-air receiver, i.e., − z rcv represents the distance from an in-air receiver to the ground surface z = 0 m) exists for the non-localized EM fields (Sun et al. 2021). Since seismic wave velocity is much lower than the speed of light, there is a high possibility that the value of k exceeds the value of ω/V em . Thus, there is a critical incident angle θ c satisfying If θ ≤ θ c , γ em will be an imaginary number and the corresponding non-localized EM fields will be radiation EM waves. If θ > θ c , then we obtain k > ω/V em and γ em is a positive real number. Thus, the depth-dependent factor exp(γ em z rcv ) represents an exponential decay suggesting the corresponding non-localized EM fields will be evanescent waves whose amplitudes decay rapidly along the normal direction of the interface.
The speed of light usually could be 5 orders of magnitude greater than the seismic velocity V sei . Consequently, the critical incident angle θ c is a very small value suggesting only the seismic waves with a nearly normal incident angle can generate radiation EM waves. We call them interfacial radiation EM waves, because they result from the seismic-to-EM conversion at an interface. Another kind of direct radiation EM waves can be directly converted from the earthquake source (Gao et al. 2014). For a receiver in non-epicentral area, both the direct and the interfacial radiation EM waves can arrive earlier than seismic P wave by several seconds or more, because they propagate with EM wave velocity. However, these radiation EM waves are so weak that they are usually invisible in the EM signals (see Fig. 3). Some 2-D numerical investigations of  Ren et al. Earth, Planets and Space (2021) 73:130 the radiation EM waves caused by motional induction effect were conducted by Gao et al. (2019).
For the non-epicentral area on the ground surface, the incident angles of seismic waves are usually large enough and evanescent EM waves are induced. Therefore, the inair EM fields shown in Figs. 5 and 6 are mainly contributed by evanescent EM waves. It should be mentioned the evanescent EM waves are also induced at the lower side of the ground surface, but they can hardly be identified unless separated from the localized EM fields.
Considering the electrokinetic effect, seismo-electromagnetic signals caused by an earthquake source have been numerically investigated in the last two decades (e.g., Pride et al. 2004;Hu and Gao 2011;Ren et al. 2012, Fig. 6 Snapshots of magnetic field B x component computed for the same times and the same area as Fig. 4 2016a; Zhang et al. 2013;Huang et al. 2015;Sun et al. 2019). The evanescent EM waves arising from the electrokinetic effect was identified for the first time by Ren et al. (2016a). Based on further numerical modellings, Ren et al. (2016b) proposed the electrokinetically generated evanescent EM waves have significant contribution to the coseismic EM signals. This idea was later adopted by Dzieran et al. (2019), who used observational data of coseismic electric signals generated by natural earthquakes to analyze the frequency-domain seismoelectric spectral ratio, i.e., the ratio of observed electric field to the seismic ground acceleration. They found the seismoelectric spectral ratio showed a decreasing characteristic for increasing frequency and this characteristic cannot be reasonably explained unless the evanescent EM waves were taken into consideration. Recently, Ren et al. (2020) further proposed evanescent EM waves can also be generated by motional induction effect or other seismo-electromagnetic coupling mechanisms. Now, this viewpoint is confirmed at least for the motional induction effect, because the generation of evanescent EM waves due to the motional induction effect can be deduced from the semi-analytical solutions (Sun et al. 2021) if one notices that the EM field expansion coefficients contain factors like exp(γ em z rcv ) (for an in-air receiver), in which γ em = k 2 − ω 2 /V 2 em ( Re{γ em } > 0 is defined when γ em has a non-zero real part). The numerical results in Figs. 5 and 6 also provide supporting evidence.

Amplitude decay characteristic
The amplitude decay characteristic of evanescent EM waves can be used to verify the validity of our computation approach. Still using the air-solid two-layer model, we compute EM fields for 16 receivers in the air, which are evenly spaced in a vertical line segment with coordinates of x = 95.76 km, y = 80.35 km and − 1.5 ≤ z ≤ 0 km. As shown in Fig. 7, all the EM components include two groups of signals showing up around 26 s and 44 s, respectively, which are the times when P and S waves arrive at the position (95.76, 80.35, 0) km. They exhibit an amplitude decay characteristic for increasing distance from the receiver to the ground surface. These EM signals actually are the evanescent EM waves, since these receivers are obviously located in a non-epicentral area. As illustrated in Appendix B, the evanescent EM waves have an amplitude decay factor f decay related to circular frequency ω , seismic incident angle θ, seismic wave velocity V sei , and the depth of an in-air receiver z rec as Further investigations are carried out to check whether the synthetic EM signals displayed in Fig. 7 have the decay characteristic described by equation (14). These signals are filtered by a frequency-domain Hanning function that can be written as where f 1 and f 2 are the cutoff frequencies. After filtering, the frequency component at f = 0.5(f 1 + f 2 ) will remain the same, while the other frequency components will be suppressed. Figure 8 shows one group of EM signals, which are induced by direct P waves and filtered by the Hanning function at 0.72 < f < 0.88 Hz. The signals in the left two plots are normalized by the maximum amplitude of all receivers' synthetic recordings. In Fig. 8(c), the maximum amplitude of each receiver is normalized by the corresponding maximum amplitude of the ground surface receiver, i.e., the one located at (95.76, 80.35, 0) km. The solid line with dot mark and the dashed line with square R filter (f ) mark denote the amplitude decay speed of the simulated E x and B y signals, respectively. The yellow line is the theoretical decay curve determined by where 8.7044 × 10 -4 is the value determined by the peak frequency 0.8 Hz, the incident angle 0.4014π (which is determined by arctan(125/40), since the ground surface receiver has an epicentral distance of 125 km and the focal depth is 40 km), and the P wave velocity of the solid half-space. Obviously, the decay curves of E x and B y are overlapped with the theoretical decay curve. A similar situation is shown by Fig. 9 for another group of EM signals, which are induced by direct S waves and whose theoretical decay curve is given by Therefore, the computed EM signals in air show an excellent agreement with the theoretically predicted amplitude decay characteristic for a single frequency and single wavenumber. This fact might be regarded as a special validation of the EM fields computed by the proposed numerical approach.

Influence of conductivity at depth
For receivers located in the shallow subsurface (e.g., at a depth of 0.1 m), which could be a common situation in field survey, the amplitude of the recorded EM fields will be affected by the conductivity of the medium in which the receivers are located, because the induction electric current σ v × B a , which acts as the source in the Maxwell's equations to generate EM fields, is evidently influenced by the conductivity σ . The amplitudes of the EM fields generally increase for a higher conductivity of the shallow subsurface. This point was already verified by Yamazaki (2012) and Gao et al. (2019). In this section, we investigate another interesting but unexamined question: how the conductivity at depth affects the waveforms of the EM fields recorded by receivers in the shallow subsurface?
The adopted model is modified from the air-solid two-layer model used in the above by enhancing the conductivity of the solid medium in the depth range of 200 ≤ z ≤ 600 m. This makes it a high conductive layer, which often exists in the Earth's crust. We use σ 2 to indicate the enhanced conductivity. Figure 10 shows the EM fields computed for a receiver located at (95.76, 80.35, 0.0001) km, i.e., in the shallow subsurface, when σ 2 is set to be 0.1, 0.3, and 0.9 s/m, respectively. The computed seismic signals are the same as those shown as yellow lines in the top three rows of Fig. 2, because the elastic properties are unchanged. The EM signals showing up around 26 and 44 s in Fig. 10 are the coseismic signals associated with direct P and S waves, respectively. All the three magnetic components (B x , B y , and B z ) show an obvious and uniform variation trend, that is, enhanced amplitude for higher σ 2 . However, the behavior of the electric field is kind of disordered. The coseismic E x signals around 26 s exhibit first a slight amplitude decrease for σ 2 changed from 0.1 S/m to 0.3 S/m and then an amplitude increase for σ 2 further changed to 0.9 S/m. For these signals, a waveform variation along with the change of σ 2 takes place and a phase reversal can be observed between the result of σ 2 = 0.1 S/m and σ 2 = 0.9 S/m. The coseismic E x signals around 44 s seemingly display a variation trend of decreasing amplitude for higher σ 2 . The signal strength of the coseismic E y component is generally enhanced for higher σ 2 , but its amplitude variation is less dramatic than those of the three magnetic components. For the coseismic E z component, neither the amplitude nor the waveform shows evident variation.
Summing up the above, in exception to the E z component, the EM fields are obviously influenced by the conductivity at depth. Such kind of influence actually results from the evanescent EM waves. For a subsurface interface of conductivity contrast, the localized EM fields are discontinuous, and thus interfacial radiation EM waves and evanescent EM waves are induced to satisfy EM a b c Fig. 9 Synthetic recordings of a electric field component E x and b magnetic field component B x , which are induced by the direct S waves and filtered by a Hanning function at 0.72 < f < 0.88 Hz. c The amplitude decay curves of the simulated EM signals (E x and B x ) in comparison with the theoretical decay curve determined by f decay = exp(14.9606 × 10 −4 z rcv ) , where z rec represents the depth of an in-air receiver fields' continuity boundary condition. The evanescent EM waves will go through amplitude decay along the normal direction of the interface and be recorded by the receiver nearby the ground surface. When the subsurface interface is close to the ground surface, i.e., the distance between them is relatively small in comparison with seismic wavelength, the evanescent EM waves will seemingly accompany seismic waves and contribute to the coseismic EM fields. Therefore, the conductivity at depth affects the coseismic EM fields through the evanescent EM waves generated at subsurface interfaces. Yamazaki (2012) considered only Rayleigh waves and also showed that the conductivity at depth can influence the magnetic field. That result can also be explained by the evanescent EM waves.

Dealing with fault rupture spreading to the ground surface
For a shallow-focus big earthquake, fault rupture spreading to the ground surface often happens and causes devastating destructions. Such case has not been considered in previous numerical investigations of seismoelectromagnetic signals (e.g., Hu and Gao 2011;Ren et al. 2012). We now adopt a model of such case and use our numerical computation approach, which combines discrete wavenumber method, peak-trough averaging method, and point source stacking method, to compute the seismo-electromagnetic signals.
The adopted model (see Fig. 11) is made up of an air layer and six solid layers consisting of s 6 , s 7 , s 8 , s 9 , s 10 , and s 11 materials, whose properties are listed in Table 2. This seven-layer model is designed by referring to CRUST 1.0 for the velocity structure at the location of (32.5°N, 130.5°E). The second and third solid layers (made up of s 7 and s 8 materials) have the same elastic properties but different electric properties. The earthquake source is a rectangular finite fault of 20 × 12.6175 km 2 penetrating the top four solid layers. Its strike, dip, and rake angles are 0°, 72°, and 65°, respectively. Distribution of the maximum slip displacement on the fault plane is displayed in Fig. 12. The fault rupture starts from the position of (0, 0, 7) km, propagates in the fault plane with a constant speed 2.7 km/s, and spreads to the ground surface. Its seismic moment is 4.65 × 10 18 N·m, which corresponds to a M6.4 earthquake. Once again, the source time function is set to be a Ricker wavelet with a peak frequency of 0.8 Hz and a time delay of 2.5 s. Referring to the data of IGRF13 (Alken et al. 2021) for the geomagnetic field at (32.5°N, 130.5°E), we set an ambient geomagnetic field with an intensity of 4.74786 × 10 -5 T, an inclination angle of 47.15°, and a declination angle of − 7.3°. The receiver is located at (76.60, 64.28, 0.0001) km. Adopting the point source stacking method, we discretize the finite fault into 100 × 60 cells. The seismic and EM signals generated by each cell are computed by jointly using the discrete wavenumber method and peak-trough averaging method. The finally obtained seismic and EM wavefields are displayed in Fig. 13. Compared with the computational result of the air-solid two-layer model (Fig. 3), the waveforms of both seismic and EM signals become more complicated, which is surely reasonable because of the multiple reflections occurring on the interfaces. The computational time of this model is about  Fig. 13 still exhibit the coseismic characteristic, which once again demonstrates the motional induction effect is likely to contribute to the coseismic EM signals recorded during natural earthquakes (Iyemori et al. 1996; Seismic vibration velocity (v x , v y and v z ), electric field (E x , E y and E z ), and magnetic induction intensity (B x , B y and B z ) computed for a receiver located at (76.60, 64.28, 0.0001) km in the seven-layer model Honkura et al. 2000Honkura et al. , 2004Matsushima et al. 2002;Ujihara et al. 2004;Tang et al. 2010).

Conclusions
In this study, a numerical computation approach is introduced on basis of the semi-analytical solutions of Sun et al. (2021) to calculate seismo-electromagnetic signals arising from the motional induction effect due to a double couple point source or a finite fault embedded in a 3-D multi-layered media. The accuracy and reliability of the proposed method is partially verified by the comparison with CGFDM for the seismic waves propagating in a multi-layered media. The existence of the evanescent EM waves resulting from the motional induction effect is identified and the influence of the conductivity at depth to the coseismic EM signals is demonstrated. Our numerical results support the viewpoint that motional induction effect is a possible explanation of the coseismic EM phenomena observed during earthquakes (Iyemori et al. 1996;Honkura et al. 2000Honkura et al. , 2004Matsushima et al. 2002;Ujihara et al. 2004;Tang et al. 2010). Based on the numerical examples shown in this paper, we conclude that the numerical computation approach proposed here provides a useful tool for the quantitative studies of seismo-electromagnetic signals. Future case studies of applying the proposed approach can help quantify how much the motional induction effect actually contribute to the coseismic EM signals.
Obviously, both the seismic and EM solutions, i.e., Eqs. (A1, A2, A3, A4, A5, and A6), are integrations with respect to wavenumber k. The order 'm' in above equations will be limited as |m|≤ 2 if a 'source-center' cylindrical coordinate system is chosen to let the point source located in z-axis (Chen 1993(Chen , 1999. The above wavenumber integrations can be regarded as summation of several inverse Hankel Transforms, because Bessel function of first kind is contained in every integrand.