**1. Introduction**

Green function for the horizontal layered seismic field is usually derived by reflectivity method, which was proposed by Fuchs and Müller [1] and extended to many other kinds [2], like reflection and transmission coefficient matrix method [3], discrete wavenumber method [4], discrete wavenumber finite element method [5], and generalized reflection transmission coefficient matrix method [6]. In the frequency domain, the Green function, derived by the reflectivity method, can be written in the Sommerfeld integral form in cylindrical coordinate system for symmetrical media.

It is well known that the numerical evaluation of Sommerfeld integral (SI) is computationally expensive due to the oscillatory and slow convergence of the integrands. To overcome this problem, several approaches have been proposed, which can be divided into two main categories: one is the approximation of the spatial domain Green functions in a closed form where no numerical integration is needed, and the other is the numerical integration of SI in conjunction with some acceleration techniques [7]. Within the first category, discrete complex image method (DCIM), which approximates the integrand of Sommerfeld integral by a series of complex exponential functions, is commonly used for the advantages of high computational efficiency, but it needs to handle the surface wave poles contributions, which not only makes the computation complicated but also brings singularity to the near region [8], and also the calculation accuracy and effective range are difficult to be accurately estimated. For the latter category, the common practice is dividing the whole Sommerfeld integral into two parts: the first part is the path to bypass the singularity; the second part, the path to infinity, is the Sommerfeld tail integral. The

**Citation:** Liu, S.; Zhou, Z.; Dai, S.; Iqbal, I.; Yang, Y. Fast Computation of Green Function for Layered Seismic Field via Discrete Complex Image Method and Double Exponential Rules. *Symmetry* **2021**, *13*, 1969. https://doi.org/10.3390/sym13101969

Academic Editors: Peng-Yeng Yin, Ray-I Chang, Youcef Gheraibia, Ming-Chin Chuang, Hua-Yi Lin and Jen-Chun Lee

Received: 23 September 2021 Accepted: 17 October 2021 Published: 19 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

finite-range integrals may readily be evaluated by the Gauss–Jacobi quadrature [9] or by the double-exponential (DE) rule [10,11]. Mosig first used DE rules to calculate Sommerfeld integrals in [11,12] which indicated its validity of suppressing endpoint singularities. The calculation of tail integral is difficult to converge due to Bessel's oscillation and slow attenuation characteristics, so this kind of method generally requires extrapolation to accelerate convergence. The WA method has shown higher levels of convergence among various extrapolation methods [13–15]. This kind of method does not need to strictly locate the position of singularity in Sommerfeld integral but only needs to ensure that the first integral path avoids all the singularities. It has good adaptability and controllable numerical accuracy; however, this depends on the number of intervals (*n*) that are chosen to evaluate the tail region. The computational time also rapidly increases as the value of (*n*) increases [16]. Given the advantages and disadvantages of these two methods, we propose a method by combining DE rules and DCIM to calculate Sommerfeld integral.

This paper first presents the Green function of a point source in a multilayer half-space in Section 3, explaining the mathematical manipulation required to obtain the solution as a Sommerfeld integral form in the frequency domain; it then describes the principle of the proposed method with DE quadrature rules and DCIM in Section 4; finally, it corroborates the correctness of the algorithm by the frequency responses obtained from the proposed approach with those where DE rules and WA partition-extrapolation are used for half-space model, and the finite element method is used for three layers model.

#### **2. Seismic Wave Equation and Green Function**

#### *2.1. Seismic Wave Equation*

The propagation of seismic wavefield in the time domain can be simplified by the following three-dimensional acoustic wave equation:

$$\nabla^2 u(x, y, z, t) = \frac{1}{v(x, y, z)^2} \frac{\partial^2 u(x, y, z, t)}{\partial t^2} + f(x, y, z, t) \tag{1}$$

where *u* (*<sup>x</sup>*, *y*, *z*, *t*), *v* (*<sup>x</sup>*, *y*, *z*, *t*), and *f*(*<sup>x</sup>*, *y*, *z*, *t*) represent displacement, velocity, and source term, respectively. *f*(*<sup>x</sup>*, *y*, *z*, *t*) = −*<sup>δ</sup>*(*<sup>x</sup>* − *xs*, *y* − *ys*, *z* − *zs*)*s*(*t*), *s*(*t*) is the wavelet, and Ricker wavelet is used in this paper; and *<sup>δ</sup>*(*x* − *xs*, *y* − *ys*, *z* − *zs*) is the Dirac function at the source point (*xs*, *ys*, *zs*).

By Fourier-transform of Equation (1), the two-dimensional acoustic wave equation in frequency domain is obtained, and therefore, Green function for the problem is defined by the following equation:

$$\nabla^2 G(\mathbf{x}, \mathbf{y}, z, \omega) + k^2 G(\mathbf{x}, \mathbf{y}, z, \omega) = F(\mathbf{x}, \mathbf{y}, z, \omega) \tag{2}$$

where, *G* denotes the Green function, *<sup>F</sup>*(*<sup>x</sup>*, *y*, *z*, *ω*) = −*<sup>δ</sup>*(*<sup>x</sup>* − *xs*, *y* − *ys*, *z* − *zs*)*S*(*ω*) is the source term in the frequency domain, *k*(*<sup>x</sup>*, *y*, *z*) = *ω*/*v* is wave number, *<sup>S</sup>*(*ω*) is Ricker wavelet in the frequency domain, and *ω* is the angular frequency.

However, the underground medium is always viscous, which leads to wave energy loss and phase change in the process of propagation. The visco-acoustic wave equation is established to better describe the propagation of the seismic waves in this viscous medium, which is the same form as Equation (2), but the complex velocity is introduced to simulate the viscous effect [17–19]. The reciprocal of complex velocity is defined as 1 *v*\$ = 1 *v* 1 − j 2*Q* [18], where *Q* is quality factor, so the complex wavenumber is set to be *k* = *ω v* 1 − j 2*Q* , j = √ −1. In this paper, the value *Q* generated by Li Qingzhong's empirical formulaisusedfornumericalsimulation[20]

$$Q = 14 \times \left(v/1000.0\right)^{2.2} \tag{3}$$

#### *2.2. Green Function in Full-Space*

The Green function of Equation (2) in homogeneous full-space can be expressed as

$$G(x, y, z, \omega) = \frac{S(\omega)e^{-ik\_1R}}{4\pi R} \tag{4}$$

where *R* = 3(*x* − *xs*)<sup>2</sup> + (*y* − *ys*)<sup>2</sup> + (*z* − *zs*)2, *<sup>S</sup>*(*ω*) is Ricker wavelet in the frequency domain, *ω* is the angular frequency, and *k*1 is the wavenumber of the medium. The above formula can be rewritten in the Sommerfeld integral form in cylindrical coordinate system:

$$\frac{S(\omega)e^{-ik\_1R}}{4\pi R} = \frac{S(\omega)}{4\pi} \int\_0^\infty \frac{m}{m\_1} e^{-m\_1|z-z\_s|} f\_0(mr)dm\tag{5}$$

where *r* = 3(*x* − *xs*)<sup>2</sup> + (*y* − *ys*)2, *m*1 = %*m*<sup>2</sup> − *k*12.

#### *2.3. Green Function in Layered Half-Space*

Consider *n* layers symmetric structure of homogeneous medium defined by interfaces located at *z*1, *z*2, ··· , *zn*−1, as shown in Figure 1. The density of each layer is *ρ*1, *ρ*2, ··· , *ρn*; and the velocity is *v*1, *v*2, ··· , *vn*. In this paper, the source is placed in the second layer.

**Figure 1.** *n* layers structure of homogeneous medium.

Each layer satisfies the acoustic Equation (2) with parameters of the Green function *G*, wavenumber *k*, velocity *v*, density *ρ*, layer thickness *h*, and quality factor *Q*, respectively. We obtain the equations as follows:

$$\nabla^2 G\_i + k\_i^{\;i} G\_i = 0, \quad i = 1, 3, \cdots, n \tag{6}$$

$$\nabla^2 G\_i + k\_i^{\;\;2} G\_{\!\!i} = -S(\omega)\delta(R - R\_0) \quad \text{i} = 2 \tag{7}$$

where *ki* = *ωvi* 1 − j 2*Q* , j =√−1.

At the interface, the pressure *Pi* = *ρiGi* as well as the gradient of the potential for the vertical direction *∂Gi ∂zi* are continuous [21]. Therefore, the following boundary conditions can be imposed on the Green function

$$\frac{\partial G\_{i}}{\partial z} = \frac{\partial G\_{i+1}}{\partial z}, \rho\_{i}G\_{i} = \rho\_{i+1}G\_{i+1}, (i = 1, 2, \cdots, n-1) \tag{8}$$

The solutions of Equations (6) and (7) can be regarded as the summation of separate up-going and down-going waves, and therefore, it can be written in the form of Sommerfeld integral in cylindrical coordinate system as follows:

$$G\_{i} = \frac{S(\omega)}{4\pi} \int\_{0}^{\infty} \left(\mathbb{C}\_{i}e^{m\_{i}z} + D\_{i}e^{-m\_{i}z}\right) l\_{0}(mr)dm, \quad i = 1, 3, \cdots, n\tag{9}$$

$$G\_{i} = \frac{S(\omega)}{4\pi} \left[ \frac{e^{-ik\_{i}R}}{R} + \int\_{0}^{\infty} \left( C\_{i}e^{m\_{i}z} + D\_{i}e^{-m\_{i}z} \right) I\_{0}(mr)dm \right], \quad i=2\tag{10}$$

where *mi* = 3*m*<sup>2</sup> − *k*2*i* , *k*2*i* = *ω*<sup>2</sup> *<sup>v</sup>*\$*i*<sup>2</sup> , *v*\$*i* = *vi* <sup>1</sup> − j 2*Qn* , *i* = 1, 2, ··· , *n*, j = √−1. In Equations (9) and (10), *emiz* and *e*<sup>−</sup>*miz* may be infinity. To maintain numerical stability, rewrite Equations (9) and (10) into

$$G\_{i} = \frac{S(\omega)}{4\pi} \int\_{0}^{\infty} \left( \mathbb{C}\_{i} e^{m\_{i}(z - z\_{i})} + D\_{i} e^{-m\_{i}(z - z\_{i-1})} \right) l\_{0}(mr) dm, \quad i = 1, 3, \cdots, n \tag{11}$$

$$G\_{l} = \frac{S(\omega)}{4\pi} \left[ \frac{e^{-ik\_{l}R}}{R} + \int\_{0}^{\infty} \left( \mathbb{C}\_{l} e^{m\_{l}(z-z\_{i})} + D\_{l} e^{-m\_{i}(z-z\_{i-1})} \right) l\_{0}(mr) dm \right], \quad i=2 \tag{12}$$

By using the boundary conditions (8), the unknowns *C*1, *C*2, *D*2, ··· , *Ci*, *Di*, ··· , *Cn*−1, *Dn*−1, *Dn* in the above formula are solved. The coefficients of the source layer are derived firstly, and other coefficients can be obtained by recursion; then, the expression of the Green function of the layered medium is obtained.

$$D\_2 = \frac{m}{m\_2} \frac{H\_2^{d\*} \left( H\_3^{u\*} e^{-m\_2|z\_2 - z\_s|} e^{-m\_2h\_2} - e^{-m\_2|z\_1 - z\_s|} \right)}{1 - H\_3^{u\*} H\_2^{d\*} e^{-2m\_2h\_2}} \tag{13}$$

$$C\_2 = -H\_3^{\mu \*} \left( \frac{m}{m\_2} e^{-m\_2|z\_2 - z\_s|} + D\_2 e^{-m\_2 h\_2} \right) \tag{14}$$

$$D\_i = \frac{D\_{i-1}m\_{i-1}e^{-m\_{i-1}h\_{i-1}}\left(1 + H\_i^{u\*}\right)}{m\_i\left(1 + H\_{i+1}^{u\*}e^{-2m\_ih\_i}\right)}\tag{15}$$

$$\mathbb{C}\_{i} = -H\_{i+1}^{\mu\_{\ast}} D\_{i} e^{-m\_{i}h\_{i}} \tag{16}$$

where  $h\_{i} = z\_{i+1} - z\_{i}$ ,  $H\_{2}^{d} = \frac{\rho\_{1}m\_{2}}{\rho\_{2}m\_{1}}$ ,  $H\_{2}^{d\*} = \frac{1 - H\_{2}^{d}}{1 + H\_{2}^{d}}$ ,  $H\_{n}^{u} = \frac{\rho\_{n}m\_{n-1}}{\rho\_{n-1}m\_{n}}$ ,  $H\_{n}^{u\*} = \frac{\rho\_{n}m\_{n-1}}{\rho\_{n-1}m\_{n}}$ ,  $H\_{n}^{u\*} = \frac{\rho\_{i+2}m\_{i+1}}{\rho\_{i+2}m\_{i+1}}$ ,  $H\_{i+2}^{u\*} = \frac{\rho\_{i+2}m\_{i+1}}{m\_{i+2}\rho\_{i+1}}$ ,  $H\_{i+2}^{u\*} = \frac{\rho\_{i+2}m\_{i+1}}{\rho\_{i+2}\rho\_{i+1}}$ ,  $i = 3, \dots, n-1$ .
