1 Introduction

Since the discovery of seismoelectric pheno-menon^{[1, 2]}, many researchers have been focused on the theoretical and experimental study of seismo-electric phenomenon^[3~15]. It should be stressed that Pride^[9] derived a set of macroscopic dynamic governing equations describing the propagation of coupled seismic and electromagnetic wave-tields in porous media, which has been widely adopted in the study of seismoelectric effect from then on.

Chen^[16~18] developed a method to calculate synthetic seismograms for layered media, decom-posing every vector using a set of orthogonal and complete basis vectors in a cylindrical coordinate system before applying the generalized reflection and transmission coefficients method. Ge and Chen^[19] developed a new efficient approach for simulating wave propagation in multi-layered media with irregular interfaces, introducing a more straightforward evaluation of the generalized reflection and transmission coefficients.

Based on the model described by the macro-scopic dynamic governing equations of Pride^[9], the above method of calculating theoretical seismo-grams^[16~18] is extended to the numerical simulation of seismoelectric wave-fields in porous media^[20]. In this new numerical emulation method, a set of scalar and vector wave-function basis is applied to decompose the concerned scalar and vector fields in the governing equations. Then two sets of governing equations corresponding to PSVTM (PSV waves generate electrical currents in the PSV particle motion plane coupling to the electro-magnetic (EM)-wavefield components with tran-sverse magnetic (TM) polarization) and SHTE (SH waves generate electrical currents in the SH particle motion plane that couple to the EM-wavefield components with transverse electric (TE) polarization) models respectively are obtained. Finally the equations are solved by the generalized reflection and transmission coefficients method, during which the straightforward evaluation of the generalized reflection and transmission coefficients developed by Ge and Chen^[19] is applied.

Compared with the previously introduced method on the numerical simulation of coupled seismic and EM waves in layered porous media^{[11, 13]}, the new technique by Ren et al. ^[20] is quite transparent for a theoretical understanding and easy on programming. Furthermore, this method provides a straightforward evaluation of the generalized reflection and transmission coefficients as well as the evaluation of the vertical component of EM waves.

However, in the evaluation of the integral of the source terms, an exponential growth factor in the integrand of integral may lead to numerical instability. The numerical calculation is stable when the source layer is thin compared with the seismic wavelength (in other words, low-frequency). However, if the source layer is thick (or at high-frequency), the numerical instability problem will emerge and prevent us from obtaining the correct results. A natural regularization approach is to introduce a fictitious thin flat source layer which includes the source point in numerical calculations. In such a numerical regularization procedure the thickness of the fictitious thin source layer is finite, which cannot always be adequate for all cases, because an appropriate thickness of fictitious layer depends on the particular situation of each case. Moreover, the introduction of fictitious source layer would result in a certain amount of extra computation. This high-frequency instability problem also exists in the calculation of theoretical seismograms. Chen^[17] introduced an analytical regularization approach which works well for any case without sacrificing any computational efficiency. In this article, we extend the analytical regularization approach of Chen^[17] to deal with the high-frequency instability problem in the numerical simulation of seismoelectric wave-fields in layered porous media. It turns out this analytical regularization approach completely solves the high-frequency instability problem. A numerical example is used to investigate the high-frequency instability problem.

2 High-Frequency Instability

The multi-layered model is assumed to have N homogeneous layers over a half-space. The j-th layer is bounded by horizontal flat interfaces, z=z^(j-1) and z=z^(j).The source placed in the s-th layer is at depth of z=z_s.

Decomposing the concerned seismoelectric fields with the scalar and vector wave-function basis, we obtain the fundamental equation given as^[20]

where j=1, 2, ···, N+1, δ_{j, s} is Kronecker delta function, y and F are the displacement-stress-EM matrix and source matrix respectively, and A is the coefficient matrix including the media properties and the electrokinetic coupling coefficient.

The general solution of equation (1) is written as

(2)

where

(3)

Θ^(j) and Λ^(j) (z) are the matrices related to the eigen-vectors and eigen-values of A^(j). The diagonal matrix Λ^(j) (z) contains exponential decay terms only. The wave-amplitude a_d^(j) and a_u^(j) are unknown vectors corresponding to down-going and up-going waves respectively. b_d(z) and b_u(z) relate to the source as follows

(4)

(5)

Through the generalized reflection and transmission coefficients and , the unknown wave amplitude vector a_d(j) and a_u(j) is determined as

(6)

and

(7)

(8)

where s_d and s_u are source terms corresponding to down-going and up-going waves respectively, and they are written as

(9)

(10)

To obtain the explicit value of wave-amplitude vector, we have to determine the source term according to equations (9) and (10). Since matrices Λ_d^(s) (ξ) and Λ_u^(s) (ξ) contain exponential decay terms only, matrices [Λ_d^(s) (ξ)]^-1 and [Λ_u^(s) (ξ)]^-1 contain exponential growth terms only, which will lead to numerical instability for high-frequency or thick source layer.

3 Analytical Regularization Approach

An efficient way to solve the high-frequency instability problem is the analytical regularization approach which creates an infinite thin fictitious source layer within the original source layer (Fig. 1).

The fictitious upper and lower interfaces are denoted by s^- and, s⁺ respectively. The correspon-ding generalized reflection and transmission coefficients at the two fictitious interfaces are given as

(11)

(12)

In obtaining the equations (11) and (12), we have used the fact that the thickness of the fictitious layer Δh→0 and

Matrix A^(j) stands for the properties in j-th layer, thus we have A^(s-)=A^(s+)=A^(s), which yields Θ^(s-)=Θ^(s+)=Θ^(s), and γ_i^(s-)=γ_i^(s+)=γ_i^(s).(13) According to the definition of Λ^(j)(z), we have

(14)

where

(15)

(16)

Substitution of equations (13)~(14) into (11) and (12) yields

(17)

(18)

which means

(19)

(20)

Under the "new" structure, the source terms s_d and s_u become

(21)

(22)

where

(23)

(24)

Since the thickness of fictitious source layer is infinitely small, the new source terms defined by equations (21) and (22) will not contain any exponential growth factors, thus the solution in frequency-wavenumber domain will be regular for any wavenumber.

According to the new generalized reflection and transmission coefticients and the new source terms, the general solutions can be rewritten as

(25)

(26)

(27)

The high-frequency instability problem is now resolved. The regularization solution given by equations (25)~(27) works well for any case.

4 Numerical Example of High-Frequency Instability

In the calculation of seismoelectric wave-fields in layered porous media, the numerical simulation is stable at low-frequencies. However, the high-frequency instability problem emerges at high-frequencies or when the source layer is thick. For each frequency, there is a critical thickness of the source layer, which separates stability from instability in numerical simulations. Now we use a three-layer-model to investigate the critical thick-ness at different frequencies.

The model is a sand layer sandwiched between two identical half-spaces (see Fig. 2). The media properties used in the numerical simulation are given in Table 1. An explosion source is positioned in the middle of the source layer whose thickness is alterable. One receiver is positioned in the top half-space 200 m above the source with horizontal offset x=10 m. The Ricker wavelet as applied as source time function. Several different center frequencies are chosen from 200 to 1000 Hz. The variation of the critical thickness of the source layer is shown in Fig. 3. The results indicated that the critical thickness decreases as the center frequency increases.

We now consider the three-layer-model with source layer thickness of 200 m and center frequency of 200 Hz. This is a thick source layer configuration. The exponential growth factor w il resuit in instability in numerical calculation. The natural regularization approach (NRA), which is introducing a fictitious thin flat source layer with finite small thickness, and the analytical regularization approach (ARA) described in this paper are respectively applied to eliminate the high-frequency instability. While using the NRA, we choose several different thickness of the fictitious source layer. It turns out we can not obtain the numerical results until the thickness is less than the critical source layer thickness (which is approximate 69 m for this configuration). The comparisons of nonzero waves in numerical results given by both approaches are displayed in Fig. 4. The amplitudes of direct P waves and accompanying electric waves are multiplied by a small factor of 0.03 to make other waves more visible. There is no significant difference between these two results. However, the computation time of ARA is less than the computation time of NRA by approximately 30%.

5 Conclusion

We discuss the high-frequency instability problem in the numerical smulation of seismoelectric wave-fields in multi-layered porous media and analyze the cause of this problem. A natural regularization approach is to introduce a fictitious thin flat source layer with a finite small thickness. In such a numerical regularization procedure, we have to find out an appropriate thickness for each case. Moreover, the introduction of a fictitious source layer would result in a certain amount of extra computation. We present an analytical regularization approach, which is creating an infinitely thin source layer within the original source layer, to eliminate the exponential growth factors in the integrand of integral of source terms. This regularization procedure does not sacrifice any computational efficiency. The results showed that the analytical regularization solution works well for any case. The technique presented in this paper may have potential applications in oll exploration^[8] and earthquake-related electromagnetic studies^[21~27].

Acknowledgments

We thank R.Robinson for reading the manuscript. This research was supported by the National R & D Special Fund for Public Welfare Industry of Earthquakes (200808069) and the National Natural Science Foundation of China (40974038, 40774028, 40821062).