*3.1. Homogeneous Model*

To verify the correctness of the generalized finite difference method for solving the elastic wave equation, a two-dimensional homogeneous medium model was set (2400 m × 2400 m), with a P-wave velocity of 4000 m/s, an S-wave velocity of 2300 m/s and a density of 2000 kg/m3. The dominant frequency of the Ricker wavelet was 20 Hz, and the source was located at the center of the model. For comparison, all sources were loaded onto the vertical component of the displacement. The generalized finite difference method was compared with both the analytical solution and the staggered grid finite difference method (SGFD). The node distribution was consistent with the rectangular grid of the staggered grid finite difference method, and the spacing of the grid (nodes) was 10 m.

To ensure the stability of the simulation, the time step was set at 0.5 ms, and the total computation time at 1.0 s. The forward modeling was performed using a 2nd-order precision 13 points GFDM, a 4th-order precision 21 nodes GFDM, and an 8th-order staggered grid finite difference method. Figure 6 shows a comparison of the 400 ms wavefield obtained by the different methods. The wavefront obtained by the GFDM is consistent with the analytical solution, which proves the correctness of the GFDM. The dispersion energy of the shear wave can be seen in Figure 6a,b, which was obtained through the 2nd-order 13 points GFDM. The results obtained using the 4th-order 21 nodes GFDM are almost the same as those obtained from the SGFD method, without any visible dispersion. This shows that improving the difference order can achieve a high calculation accuracy. In further comparing the shot records (Figure 7), the same conclusion can be drawn. In extracting the 80th trace from the shot record for comparison (Figure 8), little difference can be seen between the 4th-order GFDM results and the analytical solution, and the overall calculation accuracy is almost equivalent to that of the SGFD method. Therefore, it can be concluded that the time domain GFDM for elastic wave simulation is correct and effective.

**Figure 6.** 400 ms wavefield of homogeneous model. (**a**) 2nd-order GFDM X component; (**b**) 2nd-order GFDM Z component; (**c**) 4th-order GFDM X component; (**d**) 4th-order GFDM Z component; (**e**) SGFD X component; (**f**) SGFD Z component; (**g**) Analytical solution for the X component; (**h**) Analytically solved Z component.

**Figure 7.** Shot record of homogeneous model, (**a**) 2nd-order GFDM X component; (**b**) 2nd-order GFDM Z component; (**c**) 4th-order GFDM X component;(**d**) 4th-order GFDM Z component; (**e**) SGFD X component; (**f**) SGFD Z component;(**g**) Analytical solution for the X component; (**h**) Analytically solved Z component.

**Figure 8.** Comparison between the 80th trace of shot records: (**a**) X component; (**b**) Z component.

### *3.2. Two-Layer Model*

When SGFD methodology is used for elastic wave simulation, the grid step is determined. If the interface is located at an integral multiple of the step, the grid can accurately describe the interface; otherwise, it will cause an inaccuracy in the travel time of the reflected wave. The GFDM can directly set the nodes at the interface, therefore it can accurately describe the changes in the interface and can obtain more accurate travel time simulation results. A two-layer horizontal model (2000 m × 2000 m) was designed to show the advantages of the GFDM over the SGFD model. The P-velocity, S-velocity, and density of the upper layer of the model were 4000 m/s, 2300 m/s, and 2.4 g/cm3, respectively, and those of the lower layer were 6000 m/s, 3500 m/s, and 2.6 g/cm3, respectively. The simulation time step was set at 0.5 ms, and the total computation time at 1.0 s. The source was located at (1000 m, 0 m), and the dominant frequency of the Ricker wavelet was 20 Hz. When the model is discretized, the SGFD model adopted a fixed grid step of 10 m. Therefore, only when the interface was located at an integral multiple of 10 (for example, when the interface was located at 1000 m), can the grid points accurately describe the interface. However, when the interface is located between 1000 m and 1010 m, the results obtained by SGFD are the same, and the travel time of the reflected wave will not change. When the 4th-order 21 nodes GFDM was used, the proposed node generation algorithm was applied to discretize the model, using a node radius of 10 m for the first layer, and 12 m for the second layer. The coordinates of the nodes change along with the real interface position, so more accurate travel time information of the reflected wave can be obtained. When the interfaces were located at 1002 m, 1005 m, and 1008 m, respectively, the GFDM and SGFD results obtained are shown in Figure 9. The black and blue dotted lines represent the seismic records obtained by SGFD when the interface depth was set as 1000 m and 1010 m, respectively. The green, red, and blue solid lines correspond to seismic records obtained by GFDM when the interfaces were located at 1002 m, 1005 m, and 1008 m, respectively. Comparing the reflected wave information for the two algorithms, the SGFD can only obtain the exact record time of interfaces located at 1000 m and 1010 m—when the interface depths are 1002 m, 1005 m, and 1008 m, the record time will never change. In contrast, for the GFDM results, the record time shifts with the interface depth, which is more consistent with reality compared to the SGFD results. This shows that the GFDM can accurately describe the changes in the interface, and can obtain more accurate travel time simulation results.

**Figure 9.** Seismic records at different interface depths: (**a**) X component; (**b**) Z component.
