Stable Finite-Difference Methods for Elastic Wave Modeling with Characteristic Boundary Conditions

In this paper, a new stable finite-difference (FD) method for solving elastodynamic equations is presented and applied on the Biot and Biot/squirt (BISQ) models. This method is based on the operator splitting theory and makes use of the characteristic boundary conditions to confirm the overall stability which is demonstrated with the energy method. Through the stability analysis, it is showed that the stability conditions are more generous than that of the traditional algorithms. It allows us to use the larger time step τ in the procedures for the elastic wave field solutions. This context also provides and compares the computational results from the stable Biot and unstable BISQ models. The comparisons show that this FD method can apply a new numerical technique to detect the stability of the seismic wave propagation theories. The rigorous theoretical stability analysis with the energy method is presented and the stable/unstable performance with the numerical solutions is also revealed. The truncation errors and the detailed stability conditions of the FD methods with different characteristic boundary conditions have also been evaluated. Several applications of the constructed FD methods are presented. When the stable FD methods to the elastic wave models are applied, an initial stability test can be established. Further work is still necessary to improve the accuracy of the method.


Introduction
The propagation dynamics of seismic waves in fluid saturated porous media are of great importance for reservoir rock characterization and attract many geoscientists. The elastic waves travel through the underground material with attenuation and dispersion which is closely related to the heterogeneities of the porous continuum properties [1,2]. Attenuation is the exponential decay of the wave amplitude with distance and dispersion means the variation of wave velocity with frequency. It is commonly accepted that the wave propagation phenomenon can be described by the partial differential equations as functions of time t and space position x [1,3]. Various theories [4][5][6][7][8] focus on the dissipation mechanism of wave energy. Among the various theories, the Biot and squirt-flow mechanisms are believed to be the most important ones [6,9]. They have served as the rigorous and formal foundations to study acoustic wave propagation in saturated porous media. Numerous efforts are made to discuss the different form and numerical implementation of these two mechanisms [10]. Based on these theories, recently much attention has been paid in the literature on investigating the seismic wave propagation theories by various numerical methods, such as finite-difference (FD) and finite-element (FE) methods. When considering the simulation of wave propagation in saturated media, FD methods are well received and applied with success to numerous physical problems [11].
On the elastic wave numerical simulation or reverse-time migration, the FD method with a suitable grid is widely applied, such as the staggered grid and the rotated staggered grid. In the first-order velocity-stress hyperbolic system, the FD staggered-grid method applied with the Yee scheme is popular, but the stability of the numerical solution has been shown only in some cases [12]. Virieux provides the stability condition of the 2-order difference scheme based on a staggered grid [13]. However, when applying the FD staggered grid method on the anisotropic media, it is found that the global error significantly increases [14]. For this, Saenger et al. [15] proposed the FD rotated staggered grid method by which the accuracy is increased. Subsequently, it was discovered that this method is unstable on the boundary [16].
The construction of the algorithm usually focuses only on the accuracy. However, the stability condition is also a key factor during the elastic wave numerical simulation. The condition remands that the time step length of discrete differential should be small enough to ensure that the computation in stable. Thus far, extensive efforts have been made to overcome the limitations in calculations. Even though these conventional methods have offered great success in many applications, there can be some uncertainties and extra discrepancies owing to the lack of stable difference method and suitable boundary conditions. Thus, a stable and accurate FD method is needed to address such problems. In addition, geophysicists normally use the Fourier method to discuss the stability of the difference scheme and get some criteria, such as Von Neumann conditions [17]. This method is simple but limited to the linear differential equations with constant coefficients. For this, an energy method that can be directly applied to the nonlinear systems is considered in this work.
Some interesting combinations of difference scheme and boundary condition might overcome these limits. The idea is to use the difference scheme on the basis of the time-splitting method and to consider the adjusted characteristic boundary conditions to assure the overall stability and accuracy of the system, which is shown in Section 2. Its validity, accuracy, and applicability are systematically and theoretically investigated in Section 3. We also test the new stable finite-difference methods by applying it to the stable Biot model and unstable BISQ model with some experimental data-parameters from the paper by Yang et al. [18]. All numerical results encountered in Section 4 show that the finite difference methods and characteristic boundary conditions are suitable to simulate the wave propagation in a long time period. Therefore, for the stability analysis of different wave propagation models, the schemes are appropriate choices.

Related Work
This section presents the governing equations and the related solutions of the elastic wave modeling. The one-dimensional medium of the Biot/squirt (BISQ) elastic-wave equations can be represented as [6]: Here, v = v(x, t) and w = w(x, t) are the respective solid's and fluid's velocities at the position-time point (x, t) ∈ R × [0, T], τ is the total stress of the bulk material, P = P(x, t) is the total fluid pressure; φ ∈ (0, 1) is the porosity of the solid, M is the uniaxial modulus of the skeleton in drained conditions, α is the poroelastic coefficient, η is the viscosity of the fluid, κ is the permeability of the skeleton, F and S are the Biot-flow and characteristic squirt-flow coefficients, R = (ρ 1 ρ 22 − ρ 2 ρ 12 ) −1 , and ρ s , ρ f , and ρ a are the respective solid's, fluid's, and additional coupling densities. Subscripts x and t indicate partial derivatives (e.g., P x = ∂P ∂x , u tt = ∂ 2 u ∂t 2 ). Set U = (τ, P, v, W) T .
The linear system (1) can be written in the vector form: with coefficient matrices The previous research [19] proved that, if S > 0, there exists a symmetric positive matrix A 0 such that the matrix A 0 A 1 symmetric and A 0 A 2 is symmetric non-positive definite. Additionally, the expression of matrix A 1 shows that and the sum of its eigenvalues is zero. In this way, the matrix A 1 has two positive and negative eigenvalues. Yong [20] shows that we can diagonalize the coefficient matrix A 1 with a transformation L ∈ R 4×4 such that Here, Λ + = diag(λ 1 , λ 2 ), Λ − = diag(λ 3 , λ 4 ) and λ i (i = 1, 2, 3, 4 ) are the eigenvalues of the matrix A 1 . With the above result of the matrix A 1 , we assume that Lemma 3.3 in the paper [20] proves that there exist symmetric positive definite matrices M 1 , With these preparations, (2) can be rewritten as Setting V = LU, this can be further rewritten as where B 2 = (L T ) −1 A 0 A 2 L −1 . Additionally, the congruent transformation of the matrix does not change the numbers of positive and negative eigenvalues, and then B 2 is symmetric non-positive definite.
Corresponding to a control law [20], the characteristic boundary conditions can be prescribed as are separated for the block matrix Λ and K is a constant 4 × 4 diagonal matrix, which is defined in Section 3.

Methodology and Formulation
For the numerical integration of the system (3), the time-splitting method described in the following is considered. It means that Equation (3) are split into two separated homogeneous partial differential equations It can be assumed that a rectilinear grid with sides parallel to the coordinate axes is superimposed on {0 ≤ x ≤ 1} × {t ≥ 0} with grid spacing ∆x > 0 and ∆t > 0 in the space and time coordinate directions respectively, where J = 1 ∆x is positive integer. Let U n j = U(x j , t n ) at the grid point (x j , t n ) = (j∆x, n∆t).

First-Order Scheme
The differential Equation (5a) can be rewritten as For the interior points, it means that j ∈ [1, J − 1], the above differential equation can be discretized by , the boundary valueV n 0,+ andV n J,− can be obtained from (6) and the others can be determined by the character boundary condition (4). The next step is approximating (5b) by implicit treatment of the source term, and then U n+1 where I represents the identity matrix and then the matrix I − ∆tA 2 is invertible for the negative semidefinite matrix A 2 . The scheme (4)/(6)/(7) is a first-order finite difference schemes and the proof is given in Appendix A.1.

Second-Order Scheme
In this section, some approaches are applied to establish a second-order scheme using some thoughts guided by the research [21] and making some changes to solve (5a). Referring toV ± marked in Section 3.1, (5a) can be rewritten as where V(m) and λ m represent the m-th element in the vector V and the diagonal of the matrix Λ respectively. Obviously, λ m > 0, m = 1, 2 and λ m < 0, m = 3, 4.
In the first step of solving (8) at each time step, and, for the sake of simplicity, we set u = V(m) and obtain In the second step of the procedure, is taken to maintain the consistence of the scheme at the interior points(j ∈ [1, J − 1]) and boundary points(j = 0, J), which can be observed in calculating the truncation errors of the scheme (9)/(10)/(11) (see Appendix A.2). Adding toV n J,+ andV n 0,− calculated by the characteristic boundary conditions (4), the solutionŝ V n j of (5a) at each time step can be obtained. To maintain the difference scheme accuracy in the next step, (5b) is considered with , and the matrix M − ∆t 2B is invertible for the negative semidefinite matrixB and the positive definite matrix M.
The scheme (9)/(10)/(11) is a second-order finite difference scheme and the proof is given in Appendix A.2. Significantly, from the calculations of the truncation errors, it can be observed that the accuracy of the scheme will be destroyed if the other schemes (such as high-order Rung-Kutta difference scheme and the likes) are applied instead of the complex scheme (11) as above. In the next section, we turn to considering the stability of these schemes.

The Stability Criterion
With the energy method, we study the stability of scheme (6)/(7), and obtain the following theorem.
Proof of Theorem 1. By multiplying both sides of (6) with (LÛ n j ) T M, it can be rewritten as whereV n j = LÛ n j , V n j = LU n j . Observe that the matrices where λ m are eigenvalues of the matrix A 1 . This inequality is the CFL condition of the scheme (6)/(7). Note that M 1 Λ + and −M 1 Λ − are symmetric positive definite matrices, sum (12) over j from −∞ to ∞, and then On the other hand, by multiplying both sides of (7) with (U n+1 j ) T A 0 , it is obvious that This completes the proof.
The above theorem indicates the stability of the scheme (6)/ (7). This method is applicable to the 3D BISQ model with the same stable accuracy. Additionally, the scheme which approximates (5b) by explicit treatment of the source term remains stable with the corresponding CFL condition The proof is similar to Theorem 1. For the scheme (6)/(7) with the characteristic boundary conditions (4), Theorem 1 can be strengthened as follows.

Theorem 2.
If the time-step size ∆t, spatial increment ∆x, and the eigenvalues λ k of the matrix A 1 satisfy the CFL condition ∆t ∆x Max [|λ k |] ≤ 1 and the boundary matrix K satisfies then the solutions U n j of the scheme (4)/(6)/(7) satisfy The above theorem indicates the stability of the scheme (4)/(6)/(7) and the proof is given in Appendix B.1. This method is also applicable to the 3D BISQ model with the characteristic boundary conditions and keeps the same stable accuracy. The proof is similar to Theorem 2. Referring to the paper [20], the constant diagonal matrix K in (4) can be set as where I 2 is 2 × 2 identity matrix and κ ± are constants which follow from the inequality in Theorem 2 , such as K = 0. Additionally, for the characteristic boundary conditions the stability conditions can rewritten as follows.

Numerical Procedure and Results
In order to clarify the utility of the proposed method, the numerical simulations of the stable Biot model and the unstable BISQ model under reasonable and realistic conditions are solved and discussed (see Table 1). It is considered that the domain x ∈ [−500, 500], the spatial increment ∆x = 10, time-step size ∆t = 0.001, and the source term is a vertical point force with a Ricker wavelet with peak frequency f 0 = 35 Hz and π = 3.14. The source is set at the center of the calculation domain (x = 0). Here, length is in meters and time in seconds. T represents the calculated length in time.
The symbols and acroyms are listed with Tables A1 and A2 in Appendix C. The results using the schemes (4)/(6)/(7) and (9)/(10)/(11) are shown in Figure 1. Tables 2 and 3 show the L 1 , L 2 and L ∞ relative errors of these schemes (4)/(6)/(7) and (9)/(10)/ (11). The domain x ∈ [−10, 10], the arbitrary positive integer N, the spatial increment ∆x = 20/N, and time-step size ∆t = 0.0002∆x are considered in Tables 2 and 3.   This test has all the numerical simulations with the different model solved by the stable finite-difference methods: the stable Biot model and the unstable BISQ model. The propagation of the seismic wave for a long time period can be observed: for T = 0.05 s the quick P-wave goes ahead, in T = 0.2 s, the quick P-wave reached and was absorbed by the boundary, the slow P-wave occurs, in T = 1.5 s, the slow P-wave reached, and was absorbed by the boundary. In the front of this procedure, the seismic waves calculated by different model and finite difference method are all propagating smoothly and entirely absorbed. However, the difference occurs from T = 0.2 s. Figure 2 offers the amplitude (noted as Am) which changes as time (noted as T). By comparing (a) and (c) in Figure 1, it can be observed that the slow P-wave occurs only in the second-order scheme. The unstable phenomena only occur in the cases with the unstable BISQ model (see Figure 1b,d), in other words, the velocity occurs exponentially increasing phenomenon during the wave propagation. In the past numerical simulations, the little unstable phenomena such as those occurred in Figure 1d were often met and always attributed to the accuracy of the scheme and the boundary conditions. However, by observing Figure 1, it can be found that the instability of the numerical simulations due to the instability of the BISQ model. In addition, only for the cases in which time is large enough (such as T = 1.5 s), the unstable phenomenon which represented as the exponentially increasing amplitude as time can be significantly observed. Otherwise, the results calculated by the schemes (4)/(6)/ (7) and (9)/(10)/(11) will be stable for the stable model (see (a) and (c) in Figures 1 and 2). Even though there maybe some non-smooth solutions for the accuancy of the scheme, the amplitude also decreases to zero as time (see Figure 1a). Above all, it is pointed out that the schemes (4)/(6)/(7) and (9)/(10)/(11) are applicable for the stable model. The schemes can be used to describe the wave propagation in a long time period. The characteristic boundary condition (4) can entirely absorb the waves. It is a convenient and effective boundary condition which can be easily combined with other models for wave propagation.

Conclusions
In this work, the stable theories in mathematics are introduced into the solutions of the wave propagation theories in geology and try to reconsider the calculated results from the point of view of stability instead of accuracy. The FD method with characteristic boundary conditions that confirm the overall-stability of the schemes is proposed. The stable conditions of this scheme obtained with the energy method in the context is important. Instead of the CFL conditions by the Fourier analysis method, it can be directly generated to the nonlinear elastic wave equations. The time-splitting method and characteristic boundary conditions are utilized in this approach. The time-splitting method reduces the computational complexity and the characteristic boundary conditions which are adjusted as the interior-points schemes can confirm the overall stability. Moreover, the schemes are applied to the Biot and BISQ models, and hence the numerical simulations show that the stable Biot model or the BISQ model in stable conditions can simulate wave propagations for a long time. This fact verifies the effectiveness of the FD method as a tool to detect the stability of arbitrary elastic wave modeling. In addition, the newly introduced characteristic boundary condition can vary to fit the FD method, and then it is feasible to combine it with other schemes in future.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A.1. The Truncation Errors of the Scheme (4)/(6)/ (7) By eliminating the intermediate valuesÛ n j , the scheme (4)/(6)/(7) can be represented as By using the Taylor series expansions of U(x j , t n+1 ), U(x j+1 , t n ), and U(x j−1 , t n ) about the point (x j , t n ) the truncation error can be represented as where f n j is the exact solution value of f at the grid point (x j , t n ).
Appendix A.2. The Truncation Errors of the Scheme (9)/(10)/ (11) By eliminating the intermediate valuesV n j , the scheme (9)/(10)/(11) can be represented as For S 1,+ , by using the Taylor series expansions of in the time layer, n + 1 2 has the truncation error whereBV n j = (BV n j ) + , (BV n j ) − and (BV n j ) + , (BV n j ) − ∈ R 2×1 , and then Note that, in the numerical calculations about hyperbolic equations, normally the grid ratio is set as constant, and, for the sake of simplicity, we set ∆t = O(∆x) and obtain Similarly, we can get Then, by using the Taylor series expansions of Above all, the truncation error of the scheme (9)/(10)/(11) can be represented as .

Appendix B
This appendix shows the proof of Theorems 2 and 3.
By combining the inequalities of the interior and boundary points with the characteristic boundary conditions By combining the above inequality with the inequality of the interior points, it is deduced that Because −M 1 Λ + and M 2 Λ − are negative definite matrices, there exists a diagonal matrix K such that H 1 ≤ 0 and H 2 ≤ 0, such as K = 0. In this way, we have Similarly to the proof in Theorem 1, can be deduced and then ||U n+1 || A ≤ ||Û n || A . As above, it is proved that This complete the proof.

Appendix B.2. The Proof of Theorem 3
Theorem 3 is the same as Theorem 2 except the characteristic boundary conditions. In this way, there only offers the proof of the inequalities about the boundary points.
The others are same as the proof of Theorem 2. This completes the proof.

Appendix C
This appendix is for explaining the symbols and acronyms used in this article with Tables A1 and A2.