*2.2. Principle of Combined Compact Difference Scheme*

For the construction of the compact difference scheme, Lele extended Hermite's equation in 1992 and created the compact difference scheme as follows:

$$\begin{aligned} cf'\_{i-1} + f'\_i + cf'\_{i+1} &= a\_2 \frac{f\_{i+2} - f\_{i-2}}{4h} + a\_1 \frac{f\_{i+1} - f\_{i-1}}{2h} \\ cf''\_{i-1} + f''\_i + cf''\_{i+1} &= a\_2 \frac{f\_{i+2} - 2f\_i + f\_{i-2}}{4h^2} + a\_1 \frac{f\_{i+1} - 2f\_i + f\_{i-1}}{h^2} \end{aligned} \tag{2}$$

In Equation (2), *h* is the grid spacing, *a*, *c* are the difference coefficient matrices. *f* is the function value of node *i*. *f <sup>i</sup>* and *f <sup>i</sup>* represent the first- and second-order derivatives of node *i*, respectively; *fi*+1, *fi*+2, *fi*−1, *fi*−<sup>2</sup> represent the function values of node i successively two nodes forward and two nodes backward; *f <sup>i</sup>*+1, *<sup>f</sup> <sup>i</sup>*−<sup>1</sup> represent the first-order derivatives of node *i* successively one node forward and one node backward, respectively; *f <sup>i</sup>*+1, *f i*−1 represent the second derivative of i node one node forward and one node backward, respectively.

Following [38], the wave equation of a shear wave in a two-dimensional medium is

$$\frac{\partial^2 V\_y}{\partial t^2} = (1 + 2 \ast \gamma) v s\_0^2 \frac{\partial^2 V\_y}{\partial x^2} + v s\_0^2 \frac{\partial^2 V\_y}{\partial z^2} + \rho F\_y \tag{3}$$

where *Vy*(*x*, *z*) is the displacement component of the shear wave. *vs*0(*x*, *z*) is the shear-wave velocity in the vertical direction, and *ρFy* is the shear-wave source, that is, the concentrated force source excited in the y direction on the surface. If the input medium is isotropic, the value of the anisotropic parameter *γ* = 0.

Expanding the above shear-wave equation to represent the time of *n* + 1 and *n* − 1, we obtain

$$V\_{y(i,j)}^{n+1} = V\_{y(i,j)}^{n} + \Delta t \left(\frac{\partial V\_y}{\partial t}\right)\_{i,j}^{n} + \frac{(\Delta t)^2}{2} \left(\frac{\partial^2 V\_y}{\partial t^2}\right)\_{i,j}^{n} + \frac{(\Delta t)^3}{6} \left(\frac{\partial^3 V\_y}{\partial t^3}\right)\_{i,j}^{n} + \frac{(\Delta t)^4}{24} \left(\frac{\partial^4 V\_y}{\partial t^4}\right)\_{i,j}^{n} + O\left(\Delta t^5\right) \tag{4}$$

$$V\_{\mathcal{Y}(i,j)}^{\
u-1} = V\_{\mathcal{Y}(i,j)}^{\
u} - \Delta t \left(\frac{\partial V\_{\mathcal{Y}}}{\partial t}\right)\_{i,j}^{\
u} + \frac{(\Delta t)^2}{2} \left(\frac{\partial^2 V\_{\mathcal{Y}}}{\partial t^2}\right)\_{i,j}^{\
u} - \frac{(\Delta t)^3}{6} \left(\frac{\partial^3 V\_{\mathcal{Y}}}{\partial t^3}\right)\_{i,j}^{\
u} + \frac{(\Delta t)^4}{24} \left(\frac{\partial^4 V\_{\mathcal{Y}}}{\partial t^4}\right)\_{i,j}^{\
u} + O\left(\Delta t^5\right) \tag{5}$$

By adding these two equations, omitting the higher-order term, and substituting it into the shear-wave equation, the fourth-order accuracy difference scheme of displacement field time can be obtained as follows:

$$\begin{split} V\_{\mathcal{Y}(i,j)}^{n+1} &= 2V\_{\mathcal{Y}(i,j)}^{n} - V\_{\mathcal{Y}(i,j)}^{n-1} + (1+2\gamma)(\Delta t v)^{2} \left(\frac{\partial^{2}V\_{y}}{\partial x^{2}}\right)\_{i,j}^{n} + \left(\frac{\partial^{2}V\_{y}}{\partial z^{2}}\right)\_{i,j}^{n} \\ &+ (1+2\gamma)^{2} \frac{(\Delta t v)^{4}}{12} \left(\frac{\partial^{4}V\_{y}}{\partial x^{4}}\right)\_{i,j}^{n} + (1+2\gamma) \frac{(\Delta t v)^{4}}{6} \left(\frac{\partial^{4}V\_{y}}{\partial x^{2} \partial z^{2}}\right)\_{i,j}^{n} + \frac{(\Delta t v)^{4}}{12} \left(\frac{\partial^{4}V\_{y}}{\partial x^{4}}\right)\_{i,j}^{n} \end{split} \tag{6}$$

Equation (6) can be used to advance the time step of the shear wavefield in a twodimensional medium. The equation contains the second and fourth derivatives of the sum of displacement pairs in the *horizontal and vertical* directions, and the sum of displacement pairs in the *horizontal and vertical* directions can be differentiated by the finite difference scheme.

Chu (1998) and others constructed a combined compact difference (CCD) scheme with higher accuracy as follows:

$$\begin{cases} \begin{aligned} a\_1(f\_{i+1}' + f\_{i-1}') + f\_i' + hb\_1(f\_{i+1}'' - f\_{i-1}'') &= \frac{1}{h} \sum\_{m=1}^n c\_m (f\_{i+m} - f\_{i-m}) \\ \frac{1}{h} a\_2(f\_{i+1}' - f\_{i-1}') + f\_i'' + b\_2(f\_{i+1}'' + f\_{i-1}'') &= \frac{1}{h^2} \sum\_{m=1}^n d\_m (f\_{i+m} - 2f\_i + f\_{i-m}) \end{aligned} \end{cases} \tag{7}$$

Compared with CCD, the supercompact difference scheme needs fewer adjacent nodes to obtain the high-order accuracy approximations of the first and second derivatives at the same time step. The first and second derivatives in the above equation are coupled, which can be obtained at the same time, increasing the fidelity of the waveform.

For the CCD scheme, it is assumed that the number of longitudinal and transverse nodes of the model is *m* and *n*, and h is the size of the spatial grid. We use Equation (7) to find the spatial partial derivatives *∂*2*Vy*/*∂z*<sup>2</sup> and *∂*2*Vy*/*∂x*<sup>2</sup> in Equation (3) and then express the results as

$$AE = BV\_{\mathcal{Y}}, \ FC = V\_{\mathcal{Y}}D \tag{8}$$

where *A* and *C* are the difference coefficient matrices at the left end of the CCD Equation (7), and the sizes are 2*m* × 2*m* and 2*n* × 2*n*, respectively. *E* and *F* are the first- and second-order derivative matrices of the displacement field space to be solved, and the sizes are 2*m* × *n* and *m* × 2*n*, respectively. *B* and *D* are the difference coefficient matrices at the right end of Equation (7), and the sizes are 2*m* × *m* and *n* × 2*n*, respectively. *Vy* is the displacement field matrix, and the size is *m* × *n*. These matrices are expressed as

$$A = \begin{bmatrix} 1 & 0 & a\_1 & b\_1 h \\ 0 & 1 & a\_2/h & b\_2 \\ a\_1 & -b\_1 h & 1 & 0 & a\_1 & b\_1 h \\ -a\_2/h & b\_2 & 0 & 1 & a\_2/h & b\_2 \\ & & & & & \ddots & \ddots & \ddots & \ddots \\ & & & & & a\_1 & -b\_1 h & 1 & 0 \\ & & & & & -a\_2/h & b\_2 & 0 & 1 \end{bmatrix} \tag{9}$$

$$B = \begin{bmatrix} 0 & \sum\_{m=1}^{n} c\_m / h \\ (-2\sum\_{m=1}^{n} d\_m) / h^2 & (\sum\_{m=1}^{n} d\_m) / h^2 \\ -\sum\_{m=1}^{n} c\_m / h & 0 & \sum\_{m=1}^{n} c\_m / h \\ \sum\_{m=1}^{n} d\_m / h^2 & (-2\sum\_{m=1}^{n} d\_m) / h^2 & \left(\sum\_{m=1}^{n} d\_m \right) / h^2 \\ & \ddots & \ddots & \ddots & \ddots & \vdots \\ & & & -\sum\_{m=1}^{n} c\_m / h & 0 \\ & & & & (\sum\_{m=1}^{n} d\_m) / h^2 & (-2\sum\_{m=1}^{n} d\_m) / h^2 \end{bmatrix} \tag{10}$$

*C* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 *a*<sup>1</sup> −*a*2/*h* 0 1 −*b*1*h b*<sup>2</sup> *a*<sup>1</sup> *a*2/*h* 1 0 *b*1*h b*<sup>2</sup> 0 1 *a*<sup>1</sup> *a*2/*h b*1*h b*<sup>2</sup> ... *<sup>a</sup>*<sup>1</sup> <sup>−</sup>*a*2/*<sup>h</sup>* ... <sup>−</sup>*b*1*h b*<sup>2</sup> ... 1 0 ... 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (11) *D* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ <sup>0</sup> (−<sup>2</sup> *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> <sup>−</sup> *<sup>n</sup>* ∑ *m*=1 *cm*/*<sup>h</sup>* ( *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> *n* ∑ *m*=1 *cm*/*<sup>h</sup>* ( *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> <sup>0</sup> (−<sup>2</sup> *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> *n* ∑ *m*=1 *cm*/*<sup>h</sup>* ( *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> ... <sup>−</sup> *<sup>n</sup>* ∑ *m*=1 *cm*/*<sup>h</sup>* ( *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> ... <sup>0</sup> (−<sup>2</sup> *<sup>n</sup>* ∑ *m*=1 *dm*)/*h*<sup>2</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (12) *E* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *∂Vy ∂z* 1,1 *∂Vy ∂z* 1,2 ... *∂Vy ∂z* 1,*n*−1 *∂Vy ∂z* 1,*<sup>n</sup> ∂*2*Vy ∂z*<sup>2</sup> 1,1 *∂*2*Vy ∂z*<sup>2</sup> 1,2 ... *∂*2*Vy ∂z*<sup>2</sup> 1,*n*−1 *∂*2*Vy ∂z*<sup>2</sup> 1,*n* . . . . . . . . . . . . . . . *∂Vy ∂z m*,1 *∂Vy ∂z <sup>m</sup>*,2 ... *∂Vy ∂z m*,*n*−1 *∂Vy ∂z <sup>m</sup>*,*<sup>n</sup> ∂*2*Vy ∂z*<sup>2</sup> *m*,1 *∂*2*Vy ∂z*<sup>2</sup> *m*,2 ... *∂*2*Vy ∂z*<sup>2</sup> *m*,*n*−1 *∂*2*Vy ∂z*<sup>2</sup> *m*,*n* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (13) *F* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *∂Vy ∂x* 1,1 *∂*2*Vy ∂x*<sup>2</sup> 1,1 ... *∂Vy ∂x* 1,*n ∂*2*Vy ∂x*<sup>2</sup> 1,*n ∂Vy ∂x* 2,1 *∂*2*Vy ∂x*<sup>2</sup> 2,1 ... *∂Vy ∂x* 2,*n ∂*2*Vy ∂x*<sup>2</sup> 2,*n* . . . . . . . . . . . . . . . *∂Vy ∂x <sup>m</sup>*−1,1 *∂*2*Vy ∂x*<sup>2</sup> *m*−1,1 ... *∂Vy ∂x m*−1,*n ∂*2*Vy ∂x*<sup>2</sup> *m*−1,*n ∂Vy ∂x <sup>m</sup>*,1 *∂*2*Vy ∂x*<sup>2</sup> *m*,1 ... *∂Vy ∂x m*,*n ∂*2*Vy ∂x*<sup>2</sup> *m*,*n* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (14) *Vy* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *Vy*(1,1) *Vy*(1,2) ··· *Vy*(1,*n*−1) *Vy*(1,*n*) *Vy*(2,1) *Vy*(2,2) ··· *Vy*(2,*n*−1) *Vy*(2,*n*) . . . . . . . . . . . . . . . *Vy*(*m*−1,1) *Vy*(*m*−1,2) ··· *Vy*(*m*−1,*n*−1) *Vy*(*m*−1,*n*) *Vy*(*m*,1) *Vy*(*m*,2) ··· *Vy*(*m*,*n*−1) *Vy*(*m*,*n*) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (15)

From Equation (8), *E* = *A*−1*BVy*, the odd number behavior is *∂Vy*/*∂z*, and the even number behavior is *∂*2*Vy*/*∂z*2. Similarly, the sum can also be obtained from *∂Vy*/*∂x* and *∂*2*Vy*/*∂x*<sup>2</sup> by transposing the displacement field.

At present, this method has been applied to acoustic forward modeling [39]. However, for shear-wave reverse-time migration, it is necessary to adapt to the situation of lowtransverse wave velocity and take into account the accuracy and efficiency of calculations. Based on CCD, we introduce a combined supercompact difference scheme (CSCD), and its equation configuration is as follows [40]:

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *f <sup>i</sup>* <sup>+</sup> *<sup>a</sup>*11- *f <sup>i</sup>*+<sup>1</sup> <sup>+</sup> *<sup>f</sup> i*−1 + *a*12*h f <sup>i</sup>*+<sup>1</sup> − *f i*−1 + *a*13*h*<sup>2</sup> *f <sup>i</sup>*+<sup>1</sup> + *f i*−1 + *a*14*h*<sup>3</sup> *f* (4) *<sup>i</sup>*+<sup>1</sup> − *f* (4) *i*−1 = *<sup>a</sup>*<sup>1</sup> *<sup>h</sup>* (*fi*+<sup>1</sup> <sup>−</sup> *fi*−1) <sup>+</sup> *<sup>a</sup>*<sup>2</sup> *<sup>h</sup>* (*fi*+<sup>2</sup> <sup>−</sup> *fi*−2)... <sup>+</sup> *an <sup>h</sup>* (*fi*+*<sup>n</sup>* − *fi*−*n*) *f <sup>i</sup>* <sup>+</sup> *<sup>a</sup>*21*h*−<sup>1</sup> - *f <sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>f</sup> i*−1 <sup>+</sup> *<sup>a</sup>*22 *f <sup>i</sup>*+<sup>1</sup> + *f i*−1 + *a*23*h f <sup>i</sup>*+<sup>1</sup> − *f i*−1 + *a*24*h*<sup>2</sup> *f* (4) *<sup>i</sup>*+<sup>1</sup> + *f* (4) *i*−1 = *<sup>b</sup>*<sup>1</sup> *<sup>h</sup>*<sup>2</sup> (*fi*+<sup>1</sup> <sup>+</sup> *fi*−<sup>1</sup> <sup>−</sup> <sup>2</sup> *fi*) <sup>+</sup> *<sup>b</sup>*<sup>2</sup> *<sup>h</sup>*<sup>2</sup> (*fi*+<sup>2</sup> <sup>+</sup> *fi*−<sup>2</sup> <sup>−</sup> <sup>2</sup> *fi*)... <sup>+</sup> *bn <sup>h</sup>*<sup>2</sup> (*fi*+*<sup>n</sup>* + *fi*−*<sup>n</sup>* − <sup>2</sup> *fi*) *f <sup>i</sup>* <sup>+</sup> *<sup>a</sup>*31*h*−<sup>2</sup> - *f <sup>i</sup>*+<sup>1</sup> <sup>+</sup> *<sup>f</sup> i*−1 + *a*32*h*−<sup>1</sup> *f <sup>i</sup>*+<sup>1</sup> − *f i*−1 <sup>+</sup> *<sup>a</sup>*33 *f <sup>i</sup>*+<sup>1</sup> + *f i*−1 + *a*34*h f* (4) *<sup>i</sup>*+<sup>1</sup> − *f* (4) *i*−1 = *<sup>c</sup>*<sup>1</sup> *<sup>h</sup>*<sup>3</sup> (*fi*+<sup>1</sup> <sup>−</sup> *fi*−1) <sup>+</sup> *<sup>c</sup>*<sup>2</sup> *<sup>h</sup>*<sup>3</sup> (*fi*+<sup>2</sup> <sup>−</sup> *fi*−2)... <sup>+</sup> *cn <sup>h</sup>*<sup>3</sup> (*fi*+*<sup>n</sup>* − *fi*−*n*) *f* (4) *<sup>i</sup>*+<sup>1</sup> <sup>+</sup> *<sup>a</sup>*41*h*−<sup>3</sup> - *f <sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>f</sup> i*−1 + *a*42*h*−<sup>2</sup> *f <sup>i</sup>*+<sup>1</sup> + *f i*−1 + *a*43*h*−<sup>1</sup> *f <sup>i</sup>*+<sup>1</sup> − *f i*−1 <sup>+</sup> *<sup>a</sup>*44 *f* (4) *<sup>i</sup>*+<sup>1</sup> + *f* (4) *i*−1 = *<sup>d</sup>*<sup>1</sup> *<sup>h</sup>*<sup>4</sup> (*fi*+<sup>1</sup> <sup>+</sup> *fi*−<sup>1</sup> <sup>−</sup> <sup>2</sup> *fi*) <sup>+</sup> *<sup>d</sup>*<sup>2</sup> *<sup>h</sup>*<sup>4</sup> (*fi*+<sup>2</sup> <sup>+</sup> *fi*−<sup>2</sup> <sup>−</sup> <sup>2</sup> *fi*)... <sup>+</sup> *dn <sup>h</sup>*<sup>4</sup> (*fi*+*<sup>n</sup>* + *fi*−*<sup>n</sup>* − <sup>2</sup> *fi*) (16)

In this paper, we focused on the three-point format of Equation (16).

Using the CSCD scheme to calculate the spatial partial derivative of the equation is similar to that of CCD, except for the difference coefficient matrix. The difference coefficient matrix at the left end needs to be expanded, and its sizes are 4*m* × 4*m* and 4*n* × 4*n*, respectively. *E* and *F* are changed to the first-, second-, third-, and fourth-order derivative matrices of the displacement field space to be solved, with the sizes of 4*m* × *n* and *m* × 4*n*, respectively. *B* and *D* are the difference coefficient matrices at the right end of the equation, and the size is changed to 4*m* × *m* and *n* × 4*n*.
