*3.2. Dispersion Analysis*

In numerical simulations, if the size of the spatial grid is too large, it will cause large solution errors and produce numerical dispersion [41,42]. Therefore, the analysis of the dispersion relationship is not simply an important method to evaluate the advantages and disadvantages of numerical simulation methods, but it is also an important way to determine the size of the spatial grid [43].

If we use a simple harmonic solution *f <sup>n</sup> <sup>i</sup>* <sup>=</sup> *Aexp*[*I*(*ihk* <sup>−</sup> *kv n*Δ*t*)], we can obtain the following equations for CCD:

$$\begin{cases} \begin{aligned} &I[(a\_1+a\_1)\cos\beta+1]\beta'+[(b\_1+b\_1)\cos\beta](-1)(\beta'')^2\\ &=\sum\_{m=1}^n[(c\_m+c\_m)\cos\beta] \\ &2a\_2(I\sin\beta)\beta'+(1-2b\_2\cos\beta)(-1)(\beta'')^2=\sum\_{m=1}^n((d\_m+d\_m)\cos\beta-1) \end{aligned} \end{cases} \tag{17}$$

and the following equations for CSCD:

$$\begin{cases} \begin{aligned} \frac{1}{1} &+ \frac{187}{128} \cos(\beta) |\beta' + \frac{47}{128} \sin(\beta) \beta'' - \frac{6}{128} \cos(\beta) \beta''' - \frac{2}{768} \sin(\beta) \beta^{(4)} = \frac{315}{128} \sin(\beta) \\ -\frac{325}{64} \sin(\beta) \beta' + (\frac{69}{64} \cos(\beta) + 1) \beta'' + \frac{11}{96} \sin(\beta) \beta''' - \frac{1}{192} \cos(\beta) \beta^{(4)} &= 10(\cos(\beta) - 1) \\ \frac{315}{16} \cos(\beta) \beta' - \frac{105}{16} \sin(\beta) \beta'' + (1 + \cos(\beta)) \beta''' + \frac{1}{16} \sin(\beta) \beta^{(4)} &= -\frac{315}{16} \sin(\beta) \\ \frac{645}{4} \sin(\beta) \beta' - \frac{165}{4} \cos(\beta) \beta'' - 5 \sin(\beta) \beta''' + \frac{1}{4} (\cos(\beta) + 1) \beta^{(4)} &= -240(\cos(\beta) - 1) \end{aligned} \tag{18}$$

where *β* = *hk*, *β* = *hk* , *β* = *hk* , *β* = *hk*(3), *β*(4) = *hk*(4) Solving the equations of CCD gives,

$$\beta\_{\text{CCD}}^{\prime\prime} = \sqrt{\frac{81 - 48\cos\beta - 33\cos 2\beta}{48 + 40\cos\beta + 2\cos 2\beta}}\tag{19}$$

In this equation, *β* = *kh*, *k* is the wavenumber, *h* is the spatial gird size, and *β CCD* is the second-order modified wavenumber of the CCD scheme.

If this is substituted into the two-dimensional shear-wave equation, Equation (20) becomes

$$\beta''\_{\text{CCD}(x)} = \sqrt{\frac{81 - 48\cos\beta\_x - 33\cos 2\beta\_x}{48 + 40\cos\beta\_x + 2\cos 2\beta\_x}} \beta''\_{\text{CCD}(z)} = \sqrt{\frac{81 - 48\cos\beta\_z - 33\cos 2\beta\_z}{48 + 40\cos\beta\_z + 2\cos 2\beta\_z}}\tag{20}$$

where *β<sup>x</sup>* = *βcosθ*, *β<sup>z</sup>* = *βsinθ*, and *θ* is the angle between the wave's propagation direction and the *x*-axis. The modified wavenumber of the corresponding derivatives of CSCD can be solved by substituting the numerical solution into the Equation Group (18) obtained by the simple harmonic solution without writing the corresponding analytical expression.

The modified wavenumber solution of the finite difference scheme can then be substituted into the finite difference scheme of the shear-wave equation, yielding.

$$\begin{split} 2\cos(k\upsilon\_{\text{num}}\Delta t) &= 2 - \left(v\Delta t/h\right)^{2} \left[ \left(1 + 2\gamma\right) \left(\boldsymbol{\beta}\_{\text{x}}^{\prime\prime}\right)^{2} + \left(\boldsymbol{\beta}\_{\text{z}}^{\prime\prime}\right)^{2} \right] \\ &+ \left(v\Delta t/h\right)^{4} \left[ \left(1 + 2\gamma\right)^{2} \left(\boldsymbol{\beta}\_{\text{x}}^{\prime\prime}\right)^{4} + 2\left(1 + 2\gamma\right) \left(\boldsymbol{\beta}\_{\text{x}}^{\prime\prime}\boldsymbol{\beta}\_{\text{z}}^{\prime\prime}\right)^{2} + \left(\boldsymbol{\beta}\_{\text{z}}^{\prime\prime}\right)^{4} \right]/12 \end{split} \tag{21}$$

where *vnum* is the numerical wave velocity, *v* is the true wave velocity, Δ*t* is the time step, *<sup>v</sup>*Δ*t*/*<sup>h</sup>* is the Courant number, and *<sup>I</sup>* <sup>=</sup> √−1. Equation (21) shows that the dispersion relationship for the shear-wave equation is related to the value of spatial step, propagation speed, and time step. Based on the above dispersion relation, the corresponding *kvnum*Δ*t* can be solved after *β* is determined. The ratio of numerical wave velocity to true velocity is defined as

$$
\lambda = \frac{v\_{\text{num}}}{v} = \frac{kv\_{\text{num}}\Delta t}{kv\Delta t} = \frac{kv\_{\text{num}}\Delta t}{k\frac{v\Delta t}{h}h} = \frac{kv\_{\text{num}}\Delta t}{\frac{v\Delta t}{h}\beta} \tag{22}
$$

Ideally, if there is no numerical dispersion, then the velocity ratio *λ* is equal to one. The larger the velocity ratio, the more serious the numerical dispersion. For comparison, we calculate the dispersion relation between velocity ratio and *kh* for the central difference scheme (CFD), compact difference scheme (CD), CCD, and CSCD difference scheme. Figure 1 shows the velocity ratio curve for different *θ* for the above four methods in the isotropic (*γ* = 0) condition.

**Figure 1.** Velocity ratio curves for different numerical simulation methods. The blue curve is the traditional central difference scheme, the green curve is the traditional implicit difference scheme, the red curve is the CCD difference scheme, the blue curve is the CSCD difference scheme, and the black line is the velocity ratio constant of 1: (**a**) *θ* = 0, *α* = 0.25; (**b**) *θ* = *π*/6, *α* = 0.25; (**c**) *θ* = *π*/3, *α* = 0.25.

Figure 1a–c show the velocity ratio curves for the CCD and CSCD schemes and the other two difference schemes using different *θ*. The Courant numbers (*α* = *v*Δ*t*/*h*) are 0.25, the horizontal axis *ϕ* ∈ [0, *π*] is the product of the wavenumber and the spatial step, and the vertical axis is the velocity ratio. A velocity ratio of one indicates that the numerical wave velocity is the same as the theoretical wave velocity, which shows that the method can suppress numerical dispersion better; otherwise, it indicates that the numerical dispersion of the method is worse. It can be seen that the dispersion phenomenon of the four methods is gradually intensified with the decrease in the number of spatial sampling points. The numerical dispersion of CSCD, CCD, and CD schemes is smaller than that of the CFD scheme, as their dispersion curves are closer to one. CSCD shows the best dispersion suppression, followed by CCD.

### *3.3. The Numerical Simulation Accuracy Analysis*

To compare the numerical simulation accuracy of the models, we calculate the simulation error by simulating the two-dimensional plane harmonic initial value problem and then compare the numerical simulation accuracy of the CCD, CSCD, and CFD schemes. The initial value problem of the two-dimensional plane wave can be expressed as

$$\begin{cases} \frac{\partial^2 p\_{sh}(t, \mathbf{x}, \mathbf{z})}{\partial t^2} = v^2 \frac{\partial^2 p\_{sh}(t, \mathbf{x}, \mathbf{z})}{\partial \mathbf{x}^2} + v^2 \frac{\partial^2 p\_{sh}(t, \mathbf{x}, \mathbf{z})}{\partial z^2} \\\\ p\_{sh}(0, \mathbf{x}, \mathbf{z}) = \cos \left( -\frac{2\pi f\_0}{v} \cdot \cos \theta \cdot \mathbf{x} - \frac{2\pi f\_0}{v} \cdot \sin \theta \cdot \mathbf{z} \right) \\\\ \frac{\partial p\_{sh}(0, \mathbf{x}, \mathbf{z})}{\partial t} = -2\pi f\_0 \sin \left( -\frac{2\pi f\_0}{v} \cdot \cos \theta \cdot \mathbf{x} - \frac{2\pi f\_0}{v} \cdot \sin \theta \cdot \mathbf{z} \right) \end{cases} \tag{23}$$

where *v* is the plane wave velocity, *θ* is the angle between the propagation direction and the *x*-axis, and *f*<sup>0</sup> is the frequency of the simple harmonic plane wave.

The analytical solution for the above initial value problem is

$$p\_{sl}(t, \mathbf{x}, z) = \cos\left[2\pi f\_0 \left(t - \frac{\mathbf{x}}{v/\cos\theta} - \frac{z}{v/\sin\theta}\right)\right].\tag{24}$$

To simulate the two-dimensional shear wavefield, we specify the following: *f*<sup>0</sup> = 20 Hz and *θ* = *π*/4, the wave velocity used is 1000 m/s, the length and depth of the model are 2000 m, the length of the vertical and horizontal grids are the same, and the sampling time is 1 s. The relative error of numerical simulation is calculated under different spatial mesh lengths and time steps, and it is defined as

$$E\_{\tau}(\%) = \left\{ \frac{1}{\sum\_{i=1}^{N} \sum\_{j=1}^{N} \left[ p\_{sh}(t\_{n\prime}, x\_{i\prime}, z\_{j}) \right]^2} \sum\_{i=1}^{N} \sum\_{j=1}^{N} \left[ p\_{sh}(t\_{i\cdot j}) - p\_{sh}(t\_{n\prime}, x\_{i\cdot}, z\_{j}) \right]^2 \right\}^{\frac{1}{2}} \times 100 \tag{25}$$

where *psh<sup>n</sup>* (*i*,*j*) is the numerical solution, and *psh*- *tn*, *xi*, *zj* is an analytical solution. We then compare the relative error curves of the CCD, CSCD, CD, and CFD schemes under different spatial steps (10 m, 15 m) and fixed time steps of 1 ms (Δt = 0.001 s), as shown in Figure 2. It can be seen from Figure 2 that the relative error gradually increases with increasing spatial grid length, time step, and simulation time. When the spatial step is 15 m and the time step is 1 ms, the maximum relative error of the CCD scheme is 8%, which is much smaller when compared with the maximum relative error of the CFD scheme (39%). When the smaller spatial gird size (10 m) is adopted, the accuracy is significantly improved, and the relative error is only 0.8%. The shear-wave simulation result using the CCD scheme has high accuracy and can handle the numerical simulation of the seismic wavefield with long sampling times. The simulation accuracy of the CSCD scheme is even higher than the CCD scheme. For CSCD, when the spatial gird size is 15 m, and the time step size is 1 ms, the maximum relative error is only 0.05%, further reducing to 0.0036% when a small spatial step (10 m) is used.

**Figure 2.** Relative errors of numerical simulation for different schemes and gird size: (**a**) Δx = 10 m, Δt = 0.001 s; (**b**) Δx = 15 m, Δt = 0.001 s.
