*2.5. Boundary Condition*

In practice, the seismic wave propagates in the underground infinite medium after being excited. When the elastic wave equation is used for seismic wave numerical simulation, it is difficult to simulate the infinite medium, so it is necessary to control the size of the model within a certain range to generate an artificial boundary. If the artificial boundary is not treated with a numerical method, the elastic wave will be reflected when it propagates to the boundary, and will interfere with the wavefield inside the model. For the generalized finite difference method in the time domain, the damping attenuation boundary condition does not involve the problem of stability. By selecting appropriate attenuation parameters, a good absorption effect can be achieved.

The damping factor proposed by Cerjan [31] is as follows:

$$damp = \exp\left[-a^2(I-i)^2\right] \tag{23}$$

where *I* is the grid number of the attenuation boundary thickness, *i* = 1, 2, ··· , *I*, and *α* is the absorption coefficient. It can be seen from the formula that the selection of the absorption coefficient and the thickness of the attenuation boundary directly influence the absorption effect. The generalized finite difference method discretizes the model into nodes instead of into grids, so *I* represents the thickness of the attenuation boundary divided by the exclusion radius, *r*; and *i* represents the distance between the node and the inner boundary of the attenuation layer divided by the exclusion radius, *r*. The parameters used in this study were referenced to those used by to Li (2014) [32]. Figure 5 compares the wave field before and after adding the attenuation boundary. When there is no attenuation boundary, the energy of the wave field is covered by the boundary reflection. After adding 45 layers of attenuation boundary, the boundary reflection can be effectively absorbed.

**Figure 5.** Wave field computed by GFDM, (**a**) X component (without attenuation boundary); (**b**) Z component (without attenuation boundary); (**c**) X component (45 layers of Cerjan boundary); (**d**) Z component (45 layers of Cerjan boundary).

### *2.6. Source Term*

The source term in Equation (1) is the product of the seismic wavelet and the source spatial position function:

$$f(\mathbf{x}, z; \mathbf{t}) = s(\mathbf{t})\,\varrho(\mathbf{x}, z) \tag{24}$$

where *s*(*t*) is the seismic wavelet function, and *ϕ*(*x*, *z*) represents the impulse function. Ricker wavelet [33] and Gaussian pulse functions are used in this paper:

$$s(t) = \left[1 - 2\pi^2 f\_{\text{ff}}{}^2 (t - t\_0)^2\right] \exp\left[-\pi^2 f\_{\text{ff}}{}^2 (t - t\_0)^2\right] \tag{25}$$

where *fm* is the dominant frequency of the wavelet, and *t*<sup>0</sup> is the delay time Moreover:

$$\varphi(\mathbf{x}, z) = e^{-[(\mathbf{x} - \mathbf{x}\_0)^2 + (z - z\_0)^2]} \tag{26}$$

where (*x*0, *z*0) is the source position.

### **3. Examples**
