Consider the following 2D scalar wave equation: 1 c 2 ∂ 2 u ∂ t 2 = ∂ 2 u ∂ x 2 + ∂ 2 u ∂ z 2 + s , (1) where u and c are the wavefield and the propagation velocity, respectively, and s is the source term. With the second-order finite-difference (FD) method for the temporal discretization, the matrix form of Eq. 1 can be written as (Gao et al., 2018) follows: U n + 1 − 2 U n + U n − 1 = M U n + S ¯ n , (2) where U n + 1 and U n represent the column vectors that collect the value of the wavefields at all grid nodes at the time ( n + 1 ) Δ t and n Δ t , respectively; the column vector S ¯ n = c 2 Δ t 2 S n and S n contains entries corresponding to the source location and source time function, respectively; and the size of the update matrix M is ( N x × N z ) 2 , where N x and N z are the grid numbers of discrete points along the x- and z-directions, respectively. After applying the Fourier transform in the time domain, the left-hand side of Eq. 2 in the frequency domain can be expressed as (Gao et al., 2018) follows: ( e i ω Δ t − 2 + e − i ω Δ t ) U ˜ = − 4 ⁡ sin 2 ( ω Δ t 2 ) U ˜ , (3) where U ˜ represents the forward Fourier transform of the wavefield u . Obviously, the range of the left-hand side of Eq. 2 is from − 4 to 0. The matrix M is a semi-negative definite matrix, and its eigenvalues are non-positive real numbers (Gaffar and Jiao, 2014; Li, 2014; Gao et al., 2018; Lyu et al., 2021). Therefore, the CFL stability upper bound for Eq. 2 is obtained by requiring | ε i | ≤ 4 , where ε i is the eigenvalue of the update matrix M , and the subscript i ranges from 1 to ( N x × N z ) 2 .

The MOR method can be realized by singular value decomposition (Pereyra and Kaelin, 2008; Li, 2014) or eigenvalue decomposition (Gaffar and Jiao, 2014, 2015; Basir et al., 2015, 2018). The former method can handle a non-square matrix, while the latter method can only handle a square matrix. The update matrix M is a square matrix, and it is completely applicable to eigenvalue decomposition, whose calculation is smaller and more concise than that of the singular value decomposition. Therefore, we adopted the eigenvalue decomposition method to implement the MOR method.

To introduce the MOR method, we first performed eigenvalue decomposition on the update matrix M as follows: M = V E V − 1 , (4) where the matrix V contains the eigenvectors V i of matrix M , while E is a diagonal matrix whose entries are the eigenvalues ε i of matrix M . We used the eigenvalues ε i in matrix E to form a column vector E . Using the MOR method, Eq. 2 can be abbreviated as follows: A n + 1 − 2 A n + A n − 1 = E . ∗ A n + B n , (5) where A n and B n are the column vector of expansion coefficients for U n and S ¯ n , respectively; U n = V A n = ∑ i α i n V i and S ¯ n = V B n = ∑ i β i n V i , where the expansion coefficients α i n and β i n are the ith element of the column vector A n and B n , respectively. V i is the ith column vector of the eigenvector matrix V ; the calculation symbol “ . ∗ ” represents multiplying the corresponding elements for the column vectors on both sides.

Column vectors A n , A n − 1 , E , and B n in Eq. 5 are the input for the wavefield solving iteration, and A n + 1 is the output for the iteration. After each iteration, we can obtain the updated wavefield by the following equation: U n + 1 = V A n + 1 . (6)

The calculation of Eq. 6 can output the global wavefield values, including all the spatial grid points. If we only need to output the wavefield values of a certain trace at a fixed point, we do not need to calculate Eq. 6. For example, to output the wave value at a fixed point ( x o u t ,     z o u t ) ( x o u t and z o u t are the corner marks along the x- and z-direction in the discrete coordinate, respectively), we can use the [ ( z o u t − 1 ) × N x + x o u t ] th row vector in matrix V and multiply by A n + 1 , whose calculation amount is very small.

The input expansion coefficients B n for the source are obtained by B n = V − 1 S ¯ n . Generally, the source is located at a fixed location, that is, the column vector S ¯ n contains only one element s n that corresponds to the source; other elements are all zero. The position of the element s n determines the location of the source, and the value of the element s n represents the value of the source time function at time n Δ t . Source column vectors S ¯ n + 1 and S ¯ n at adjacent moments only have one different element with a scalar factor ratio   α (e.g., S ¯ n + 1 = α S ¯ n ), which is determined by the amplitudes of the source time function at different times. Therefore, there is no need to perform the matrix calculation B n = V − 1 S ¯ n for each iteration. We only need to know the ratio α of the source time function at adjacent moments, and the expansion coefficients B n + 1 can be obtained by B n + 1 = α B n .

We can see that only the eigenvalues   E of the updated matrix and the expansion coefficients A n , A n − 1 , and B n (wavefield and source) in Eq. 5 are used to do the iteration of wave propagation, which is just an iteration of the column vectors and greatly simplifies numerical simulation compared with Eq. 2.

We used the fourth-order FD method for the spatial discretization of M in Eq. 2. Using the velocity c = 4,000 m/s and the spatial grid interval Δx = Δz = h = 10 m, we can obtain the maximum time step Δ t max = 1.530 ms, according to the CFL stability condition for high-order FD schemes (Liu and Sen, 2009). For a time step Δt over Δ t max , the discrete update matrix M will contain unstable eigenvalues (Li et al., 2014; Gao et al., 2018; Lyu et al., 2021). As shown in Figure 1, the eigenvalues of the discrete update matrix for Δ t = 1   ms are all distributed in the range of 0 to −4, which means that all the eigenvalues are stable (Gao et al., 2018); in contrast, some eigenvalues of matrix M for the time step larger than Δ t max (e.g., Δt = 2, 3, 4, 5, 6, 7, 8, and 9 ms) are distributed outside of the range of 0 to −4 (the red zones showed in Figure 1). These eigenvalues, whose absolute values violate the basic requirement of | ε i | ≤ 4 , cause unstable phenomena when Δ t > Δ t max . Table 1 shows the number of stable eigenvalues for different time steps. We can see that after the time step exceeds Δ t max , the number of stable eigenvalues decreases sharply as the time step increases (e.g., when Δt = 9 ms, only 965 stable eigenvalues exist, which is 2.39% of the total numbers 40401). To release the CFL stability upper bound on the time step, we introduced two kinds of operations for those unstable eigenvalues of matrix M : the eigenvalue abandonment algorithm and the eigenvalue perturbation algorithm.

We performed numerical experiments on a homogenous square model by the finite-difference time-domain method. The wave velocity is c = 4,000 m/s. The spatial grid interval is Δx = Δz = 10 m, and the grid point number is 201 × 201. The source is a Ricker wavelet with a dominant frequency of 20 Hz, which is located at x = 1.0 km and z = 1.0 km. According to the CFL stability condition, the maximum time step for the second-order FD method in temporal discretization and fourth-order FD method in spatial discretization is Δ t max = 1.530 ms. We tested several large time steps (Δt = 2, 3, 4, 5, 6, 7, 8, and 9 ms) by Eqs. 9–12, where all time steps are beyond the CFL stability upper bound. We used the workflow in Figure 3 to apply the time-dispersion transform method. To examine the results obtained by the different time steps, we performed numerical simulations using a short time step Δt = 1 ms, which is also applied by the time-dispersion transform method. The simulated waveforms can be regarded as theoretical references to examine the accuracy using larger time steps.

Figure 4 shows the waveforms recorded at x = 700 m and z = 700 m. Although the time steps exceed the CFL stability condition, no instability arises after applying the eigenvalue abandonment (shown in Figure 4A) and the eigenvalue perturbation (shown in Figure 4B). It proves that the combination of the MOR, eigenvalue abandonment (and eigenvalue perturbation) algorithm, and time-dispersion transform methods can extend the CFL stability upper bound successfully.

The time-dispersion error using Δt = 2, 3, 4, 5, and 6 ms is invisible, as shown in Figures 4A,B. It indicates that integration of the MOR, eigenvalue abandonment (or eigenvalue perturbation) algorithm, and the time-dispersion transform method can provide highly accurate simulation results even when a much larger time step size beyond the CFL stability upper bound is used.

According to the Nyquist sampling theorem Δ t < 1 / ( 2 f max ) (Gaffar and Jiao, 2014; Gao et al., 2018), the time step size cannot be larger than 6.7 ms for the Ricker wavelet with a dominant frequency of 20 Hz, whose maximum frequency is about 75 Hz. Therefore, using Δt = 7 ms has exceeded the Nyquist sampling. Interestingly, starting from Δt = 7 ms, the results of the two processing methods are different. Although 7 ms has exceeded the Nyquist sample for the eigenvalue abandonment operation, the results still seem acceptable except for some slight error. We need to associate the eigenvalue abandonment with the spatial filtering method (Gao et al., 2019). Different eigenvalues of the updated matrix correspond to different wavenumbers: the unstable eigenvalues correspond to the wavefield that distributes in the high-wavenumber region, and stable eigenvalues correspond to the low-wavenumber region (Li, 2014). Abandoning the unstable eigenvalues is equivalent to filtering out the unstable wavefield components with a low-pass filter. The aliasing effect caused by insufficient sampling points (i.e., using a time step that exceeds the Nyquist sampling upper bound) is a high-frequency oscillation, which would be filtered out by a low-pass filtering operation.

In contrast, the eigenvalue perturbation operation is implemented by perturbing the unstable eigenvalues into stable eigenvalues, which is a normalized operation, rather than a low-pass filter. The high-frequency oscillation caused by insufficient sampling points will be retained. It can explain why the high-frequency oscillation exists resulting from the eigenvalue perturbation algorithm using Δt = 7 ms in Figure 4B.

Next, using the idea of spatial filtering, it is easy to explain the inaccuracy for the result obtained using Δt = 9 ms in Figure 3A. As the time step increases, the number of stable eigenvalues decreases dramatically (as shown in Table 1). It is equivalent to a sharp decrease in the threshold of the low-pass filter for spatial filtering. The effective seismic wavefield mainly distributes in a certain bandwidth of low wavenumber. The corresponding low-pass filter would filter out the effective wavefield distributes in the low-wavenumber regions for an excessively large time step. As a result, the wavefield information is incomplete. Even if the time-dispersion transform method is used, the result would still have errors.

We further analyzed the amplitude errors between the waveforms obtained by different time steps and the theoretical waveform (as shown in Figure 5). The time ranges from 3 to 3.1 s, which is an intermediate time period of the waveforms in Figures 4, 5C,D, are logarithmic displays for the absolute values of the amplitude errors in Figures 5A,B, respectively. The amplitude errors increase with the increasing time step, but the errors are acceptable when Δt ≤ 6 ms. For example, using Δt = 2 ms, the values of amplitude errors are around 0.001 (red lines in Figures 5C,D), while the amplitude values of the theoretical waveform range from −3.340 to 4.011, and the error of 0.001 is negligible relative to the overall amplitude values; using Δt = 6 ms, the maximum error is about 0.1 (grey lines in Figures 5C,D), which is only 2.5% of the maximum amplitude value 4.011.

Based on the analysis mentioned before, for a homogenous model, in association with the MOR method and the time-dispersion transform method, both the eigenvalue perturbation algorithm and the eigenvalue abandonment algorithm can release the time step toward the upper bound of the Nyquist sampling and still ensure the accuracy of the numerical simulation. For the eigenvalue abandonment operation, although no high-frequency oscillations appear after the time step size exceeds the Nyquist sampling, there is a risk of filtering out effective wavefield that distributes in the low-wavenumber region if the time step size is too large.

We verified the feasibility of the proposed methods by a heterogenous medium, part of the Marmousi model, as shown in Figure 6. In this model, the velocity contrast is strong, the multiple waves are significant, and the velocity range is from 1,467 to 5,928 m/s. The spatial grid interval is Δx = Δz = 10 m, and the grid point number is 121 × 201. The source is a Ricker wavelet with a dominant frequency of 15 Hz, located at x = 1.0 km and z = 0.6 km. The maximum time step size for the second-order FD method in temporal discretization and the fourth-order FD method in spatial discretization schemes are Δ t max = 1.032 ms. We tested several large time steps (Δt = 2, 3, 4, 5, 6, 7, 8, and 9 ms), and all these time steps are beyond the CFL stability upper bound. We used the workflow in Figure 3 to apply the time-dispersion transform method. Similar to what has been done in the heterogenous model, we performed numerical simulations using a small time step Δt = 1 ms, whose results can be regarded as theoretical references to examine the accuracy using larger time steps.

Figure 7 shows the waveforms around 3 s. The residual errors are invisible even for the eigenvalue perturbation algorithm, even for Δt = 7 ms (shown in Figure 7B). Starting from Δt = 8 ms, the simulation results become inaccurate due to the time step size has exceeded the upper bound of the Nyquist sampling. It demonstrates that the combination of the MOR, the eigenvalue perturbation, and the time-dispersion transform method can still release the time step to the upper bound of the Nyquist sampling even for the heterogenous media with strong velocity contrast.

For the eigenvalue abandonment algorithm, when Δt = 6 ms, a relatively obvious error appears. Simultaneously, the time step size has not reached the upper bound of the Nyquist sampling yet (shown in Figure 7A). This problem still needs to be explained from the perspective of the spatial filtering method. In heterogenous media, instability phenomenon would appear in the high-velocity region according to the CFL stability condition with the time step size increasing. Therefore, the threshold of the low-pass filter is determined by the high-velocity region. However, for the wavefield with a given bandwidth, the wavenumber range in the low-velocity region is wider than that in the high-velocity region (Gao et al., 2019). With the increasing time step size, the low-pass filter of the spatial filtering with decreasing threshold will filter out the effective wavefield in the low-velocity region, which would result in incomplete wavefield components. For the eigenvalue abandonment algorithm, when the time step size is too large (e.g., Δt ≥ 6 ms in this numerical experiment), the unstable eigenvalues are caused by the velocity in the high-velocity region, and this part of eigenvalues corresponds to the wavefield in the high-wavenumber region. Abandoning these unstable eigenvalues is equivalent to filtering out the wavefield of the corresponding wavenumber range, which would filter out the effective wavefield of the low-velocity region.

We further analyzed the amplitude errors between the waveforms obtained by different time steps and the theoretical waveform (as shown in Figure 8). The time ranges from 3 to 3.1 s, which is an intermediate time period of the waveforms in Figure 7. The amplitude values of the theoretical waveform range from −3.345 to 3.081. Using Δt = 6 ms, the maximum error of the eigenvalue abandonment algorithm is 0.204 (grey lines in Figures 8A,C), which is 6.1% of the maximum absolute amplitude value 3.345; while the maximum error of the eigenvalue abandonment algorithm is 0.059 (grey lines in Figures 8B,D), which is only 1.8% of the maximum amplitude value 3.345. Therefore, it is better to perturb this part of the eigenvalues to be stable than to abandon them directly, and this can retain the effective component of the wavefield. The eigenvalue perturbation algorithm is a better choice than the eigenvalue abandonment algorithm. The former is more suitable for the simulation of strongly heterogenous models. The time step can be released beyond the CFL stability upper bound and even toward the upper bound of the Nyquist sampling.