Abstract

A stable and accurate finite-difference discretization of first-order elastic wave equations is derived in this work. To simplify the origin and proof of the formulas, a symmetric matrix form (SMF) for elastic wave equations is presented. The curve domain is discretized using summation-by-parts (SBP) operators, and the boundary conditions are weakly enforced using the simultaneous-approximation-term (SAT) technique, which gave rise to a provably stable high-order SBP-SAT method via the energy method. In addition, SMF can be extended to wave equations of different types (SH wave and P-SV wave) and dimensions, which can simplify the boundary derivation process and improve its applicability. Application of this approximation can divide the domain into a multiblock context for calculation, and the interface boundary conditions of blocks can also be used to simulate cracks and other structures. Several numerical simulation examples, including actual elevation within the area of Lushan, China, are presented, which verifies the viability of the framework present in this paper. The applicability of simulating elastic wave propagation and the application potential in the seismic numerical simulation of this method are also revealed.

1. Introduction

Partial differential equations (PDEs) are widely used in mathematics, physics, and other fields. However, in practical engineering, it is difficult to obtain the corresponding analytical results. Hence, a number of numerical methods have been developed for computational modeling of elastic wave propagation. For example, finite-difference methods (FDMs) [13], finite element methods (FEM) [4, 5], spectral element methods (SEM) [6], boundary element methods (BEM) [79], finite volume methods (FVM) [10], and pseudospectral methods [11, 12]. However, each method has its own merits and drawbacks.

FDMs are commonly used and extremely efficient for wave propagation in seismic engineering applications due to their efficiency and ability to simulate complex layered media, among others. However, it is difficult to manage complex geometries, particularly when using a staggered grid method [13]. An additional problem of this method is that it is unsatisfactory when simulating cracked media. Generally, when the traditional FDM is used, the method to simulate the crack is to fill the crack with a low velocity medium such as water. Meanwhile, grids close to a low velocity medium must be refined for stability, which requires a high volume of calculations and difficult automatic modeling. Additionally, selected FDM methods present difficulties in terms of managing long-time stability, both in terms of spatial- and time-related aspects. These drawbacks give rise to significant limitations for the simulation of large-scale, long-time, and undulating media. Therefore, a numerical method with time and spatial stability is proposed in this work.

In the context of numerical simulation, accuracy will have a significant effect on results. Compared to first- and second-order methods, higher-order FDM (HOFDM) can effectively reduce the degree of freedom required for simulating wave propagation [14]. However, HOFDMs also give rise to additional limitations and difficulties, such as discrete approximation of boundary conditions, and the establishment and simulation of complex structures such as cracks. Furthermore, a larger amount of computation and longer computation time are generally required, leading to computation limitations.

In large-scale simulation cases particularly, the computational load and memory utilization of traditional methods will face significant challenges. Multiblock simulation methods proposed for composite domains can save processing memory by partitioning the region represented in a large-scale model, thereby greatly improving computing efficiency by enabling parallel calculations between blocks.

A relatively appropriate HOFDM that combines initial boundary value conditions involves summation-by-parts (SBP) operators [1519]. This approach is employed for approximating the derivatives in a spatial domain, together with the weak enforcement of boundaries and interface conditions, using either the simultaneous-approximation-term (SAT) method [20, 21] or the projection method [2224]. The SBP-SAT method can develop strictly stable proof using the energy method, leading to robustness in spatial and time domains. Some review papers of SBP operators can be found in [25, 26] and examples of SBP-SAT can be found in [2730].

The essential characteristics of SBP operators are that they mimic the integration-by-parts (IBP) property so that it is possible to mimic the energy dissipation of a continuous problem. The basic theory of first-order SBP operators was proposed by Kreiss and Scherer [19] in 1974, whereby the energy method is applied for stability proof, leading to the HOFDM known as strict stability [31, 32].

The SBP operators introduced above were derived in a uniform grid with central-difference stencils in the interior of the area and using one-sided difference stencils near the boundaries, thereby maintaining the properties of SBP operators. These operators can be referred to as classical SBP operators. Other types of SBP operators are generalized SBP [25, 3335], multidimensional SBP [36, 37], upwind SBP [38, 39], and staggered and upwind SBP operators [40]. Generalized SBP operators have one or more characteristics of nonrepeating interior operators, nonuniform nodal distribution, and exclude one or both boundary nodes. Multidimensional SBP operators remove the limitations of classical SBP operators for applying the tensor-product formulation to multidimensional conditions because most classical SBP operators are one-dimensional. Upwind SBP operators, combining flux-splitting techniques for hyperbolic systems were proposed via noncentral-difference stencil in the interior, which introduced artificial dissipation on a nonstaggered or staggered grid.

In seismology, and in a number of fields such as acoustics, electromagnetics, and oceanography, a commonly encountered difficulty is the treatment of unbounded domains. The local high-order absorbing boundary condition (ABC) [4143] and perfectly matched layer (PML) [4447] are the most commonly used techniques for trimming these domains, and “absorbing” outgoing waves as much as possible, in order to simulate elastic wave motion in finite computation domains. ABC is a boundary condition defined on an artificial boundary in order to make small, spurious reflections occur when elastic waves propagate at the boundary. Absorbing layers were proposed by modifying the underlying equations, in order to rapidly dampen the wave in a layer. Due to the difficulties associated with the transformation between time and frequency domains, energy estimation was not applicable to the PML within the time domain [48]. A local high-order ABC, as well as a nonreflecting boundary condition (BC), has been widely used for their feasibility, intuitive nature, and generalizability when combined with the SAT method, enabling well-proven approximates [4958]. We employed the nonreflecting BC since it provided a simple formulation of the elastic wave problem. As an extension, a method combining PML with the SBP-SAT (summation by parts + simultaneous approximation terms) methodology of the symmetric matrix form (SMF) of elastic wave equations was developed in a general curvilinear geometry and is presented in this work. Nonreflecting BCs derived in the Laplace space can be observed in [59]. Additionally, a well-proposed SBP-SAT methodology was used on grids with nonconforming interfaces [6064].

The aim of this work was to establish a generic representation of the elastic wave equation using the SMF. As a result, the SMF framework was able to formulate SH and P-SV wave equations in a similar manner, thereby simplifying this calculation, while also facilitating the transformation of the extension of dimensions (two-dimensional [2D] and three-dimensional [3D]). Meanwhile, using the SMF, a general continuous stability analysis is presented, combined with BCs in a general hyperbolic system. An equivalent form can be derived in a curve domain, similar to the form in a cuboid domain, by using curvilinear transformation of the elastic wave equations. Thereafter, a generic stable discretization of the SMF on curvilinear coordinates was developed using SBP-SAT, where it was observed that the SMF is able to provide an equivalent type of the cuboid domain. In this instance, the results of similar stability analyses greatly simplified the calculation and proof process. Hence, this framework is applicable to a wide range of elastic wave equations.

The remainder of this paper is presented as follows: in Section 2, the SMF of elastic wave equations is established. The well-posed BCs of the SMF of elastic wave equations are proposed using the flux-splitting form provided in Section 3. Section 4 presents the stability analysis using the continuous energy method. A semidiscrete and applied SBP-SAT methodology is presented, and stability analysis within a discrete norm is discussed in Section 5. Numerical studies are performed in Section 6. Subsequently, conclusions are drawn and future research is presented in Section 7.

2. SMF of Elastic Wave Equations

2.1. Formulation of SH Wave

Let the Cartesian coordinates of a two-dimensional spatial domain defined by and t denotes the time variable. Consider the first-order elastic wave equation of the SH form:where is the velocity in the z direction of the SH wave, defined by , and are the stresses; is the density of the elastic medium; and is the shear modulus. Equation (1) can then be revised as a symmetric system that has a formwhere

2.2. Formulation of P-SV Wave

The P-SV from the elastic wave equation can be written aswhere are the particle velocities, are the stresses, is density, and is the stiffness matrix. In the context of the isotropic case, the stiffness matrix is described by two independent elastic coefficients, the Lame parameters and . Note that C is a symmetric positive definite:

It will be convenient to use wave velocities and as medium properties:

Equation (4) can be formulated as a symmetric system aswhere ,

The matrices in (8) are symmetric, and is positive definite. Equation (7) also can be written as a symmetric first-order hyperbolic system, as a generic system (2):

In order to transform (7) to (10), we shall ansatz , where

P is an orthonormal matrix that satisfies , denoting

Thus , and premultiplying M from Equation (10) yieldsas

Here, we transform (14) into

Matrix M satisfying the above conditions is

2.3. Curvilinear Coordinates and Coordinate Transformation

Regular grids have difficulty in managing numerical simulations for geometrically complex models. Consider the elastic wave equation in a domain with a curved boundary shape, where a transform from the curvilinear domain to the cuboid domain was effected (Figure 1). Here, we introduced mapping defined by

The Jacobian determinant is

Therefore, the derivatives in each direction can be written as

Using the relations, it follows that

Thus, we obtain the same formula as (2) and (10):where are symmetric matrices since and are symmetric as well.

3. Well-Posed Boundary Conditions

A generic representation of the elastic wave equation, referred to as the SMF, is presented. Using this framework, a specific class of the elastic wave equation can be analyzed and approximated. Meanwhile, the implementation of BCs is also simplified. Wave propagation is often extremely large, even in unbounded spatial domains. Considerable physical domains must be replaced by suitable-sized computational domains, by introducing artificial BCs to trim these domains. Moreover, a free-surface BC is necessary for simulating elastic wave propagation at the surface or crack structure. The traditional approach for simulating different BCs is to employ the properties of waves at different boundaries, e.g., the traction perpendicular to the surface of the free boundary is zero. Different BCs have different mathematical formulas, which gives rise to difficulties in the derivation and unity of the elastic wave equations framework. The SMF is not simply a generic representation of the elastic wave equation but also creates a generic system of different BCs. To derive a nonreflecting BC and the free-surface BC, a well-proven characteristic type of ABC was considered, where the number of variables propagated into the domain was equal to the number of variables propagated out the domain [65, 66]. Here we denote four boundaries in a domain: west, east, south, and north (Figure 1). The free-surface boundary condition at the north boundary was considered. A nonreflecting BC was proposed at the west, east, and south boundaries. More details of characteristic BCs are provided in [67, 68]. In this paper, we first introduce the splitting of the coefficient matrix:

Then, having

Using (23), the characteristic variables associated with propagating in and out of the domain at the boundary are

Well-posed BCs posit that only incoming variables of the correct quantity must be prescribed on the boundary. Thus, the well-posed problem can be specified by the incoming characteristic variables composed of T, which iswhere R serves as a reflection coefficient, i.e.,. Aimed at different BCs, R = −1 represents a free-surface BC and a nonreflecting BC, as defined by R = 0.

4. Continuous Energy Analysis

In this section, an energy analysis is presented for the SMF of the elastic wave equation (10), alongside the BC in (25). The energy method was based on constructing a norm for the given problem that did not grow over time. To simplify notation in the forthcoming energy analysis, a number of definitions are necessary. The result at the grid point was stacked as a vector, where and . The inner product was given by and in the domain . Thus, the corresponding norm was expressed as . Here, we again included the SMF of elastic wave equations in a bounded domain for the sake of convenience:

Multiplying (26) by , setting data to zero and adding its transpose, then integrating over the domain yield

Equation (27) is equivalent towhere is the time derivative, and (28) gives

A bounded energy for (28) required , as well as . In this instance, we show that the BC terms (29) were nonpositive definite. As BCs were the same except at the north boundary, we only included the treatment of the south and north boundaries:

Equations (30) and (31) were nonpositive so that (28) satisfied the energy estimate. Using the SMF, the energy analysis could be simplified to a large extent as a generic system.

5. Semidiscrete Approximations

In this section, an energy-stable, semidiscrete approximation was derived. We used SBP finite-difference operators to approximate the spatial domain and imposed BCs using the SAT method, on an unstaggered grid.

5.1. Spatial Discrete Operators

To begin with, consider the uniform discretization of the domain, , with N + 1 grid points and step h > 0:

The solution is stacked as a vector , and its discrete form is . An SBP operator approximating has the form

The matrix D satisfies integration by parts if the following properties hold:where . In this paper, 6th order SBP operators are used, as mentioned in [15].

5.2. SBP-SAT Methodology of SMF of Elastic Wave Equations

In this section, an SBP-SAT methodology of SMF of elastic wave equations is proposed and its numerical stability is proved. A computation domain was discretized using grids. The numerical solution of SH wave equations on the grids is stacked as vectors:

For the P-SV wave equations, (35) then take the form

The SBP method in higher space dimensions is derived using the Kronecker product:where B is an matrix and C is a matrix so that (37) is a matrix. (37) fulfills the rules of and . Define be the identity matrix for the SH wave equation and identity matrix for the P-SV wave equation, be the , and be the matrix; SBP operators in 2D are

The coefficient matrices are determined in 2D yields:

Moreover, the inner product is defined as follows:

To combine with boundary conditions, we set

Here, the SBP-SAT discretization of (21) with the BCs (24) and (25) can be written aswhere

For stability analysis, multiplying (42) by and adding the transpose lead towhere the boundary conditions are

Similar to (30) and (31), we have

Since (46) are nonpositive definite matrices, verifies the robustness.

For large scale, as well as complex structural spatial domains, a computation block is not a preferable choice for numerical simulation. We divided the computation space into several blocks, with grids conforming at the block interface when interfaces were continuous, which had a form matching that in (47). Meanwhile, the multiblock SBP-SAT methodology of elastic wave equations was able to manage discontinuous structures such as cracks by using the interface definition. The form of two blocks is shown in Figure 2.

when the interface conditions were applied for simulating cracks that had

6. Numerical Experiments

6.1. Multiblock Condition of P-SV Wave

In this section, we consider a multiblock condition that can divide the computation domain into appropriate scale blocks. Each block was able to transform data through an adjacent boundary, which had the ability to effect computations in parallel. To demonstrate the potential of this multiblock treatment, we imposed nonreflecting BCs at boundaries, and continuous interface conditions for blocks. The properties in blocks were  = 1800 m/s,  = 1000 m/s, and . The computation domain was bounded by , which involved nine blocks with , as shown in Figure 3. The time step is set to , and the initial data are

Using the SBP-SAT methodology for the SMF of elastic wave equations, we obtained amplitude diagrams for different times. The x- and y-direction wave fields showed different waveforms in the P-SV mode since the initial data had been applied in the x direction. P-wave propagated from block five to other blocks at time = 0.3 s (Figure 4). This indicated that the P-SV wave field motivated in block five showed little change when propagating to other blocks.

P-wave propagated to the nonreflection boundary at time = 0.9 s (Figure 5), and SV-wave propagated to the nonreflection boundary at time = 1.5 s (Figure 6). This indicated the feasibility of unbounded domain simulation.

The treatment of BCs is always more complex when applying the higher-order difference method. Generally, the accuracy of BCs is lower than that for the inner medium. Therefore, the stability problem often occurs at the boundary position. This problem becomes more prominent when multiple computational domains are combined. When doing so, we divided the domain into nine subblocks and 36 boundaries, which posed a significant challenge to the stability and impact of the numerical simulation. Stability is demonstrated in Section 5, as part of this numerical experiment. Meanwhile, when the elastic wave propagated to the inner boundary of the medium, there was no disturbance field; that is, the block form had little influence on the wave field, which indicated that the multiblock simulation of the wave equation was feasible. The SBP-SAT methodology applied in this work clearly had a positive promoting effect on numerical simulations, which can significantly expand the simulation scale.

To investigate the connection between the multiblock and propagation of the elastic wave, we next present the amplitude of velocity in the x direction at time = 0.6 s in block two independently (Figure 7). We verified that continuous interface conditions had little impact on the overall wave field when waves propagated to the interface boundaries. Meanwhile, the values of the boundary nodes of the two adjacent computational domains coincided with one another (Figure 8). It is shown that a multiblock domain can enlarge the simulation scale and significantly improve computational efficiency under reasonable settings.

6.2. Crack Structure for SH Wave

When the inner computational domain was considered as a whole, the boundary types in the subcomputational region could also be diversified, which provided convenience for simulating complex structures. For example, the traditional finite-difference method is required to perform complex operations or approximations in order to realize a discontinuous medium; for example, filling intermittent parts with low-speed media. In order to visually demonstrate the propagation process of elastic waves in cracked media, SH waves were utilized since waveform transformation was not involved. The properties of the medium were  = 100 m/s, , and . The computational domain was bounded by , which involved nine blocks with . The time step was set to , and the initial data were similar to that in (49).

When an SH wave propagated to the crack structure, it generated reflected and diffracted waves, which caused continuous disturbance of the wave field in the medium. The superposition and mutual influence of the wave field had a significant impact on wave propagation characteristics. Meanwhile, the amplitude of diffracted waves decreased in an obvious manner (Figure 9). This helped us to explain that in the case of real earthquakes, seismic wave attenuation factors must be considered based on the impact of cracks. Moreover, the method proposed in this paper is more convenient for realizing a crack structure that is able to avoid the amplification of errors by low-speed media filling.

6.3. Curve Crack Structure in Curvilinear Domain

In this numerical experiment, we considered curve crack structure in a curvilinear domain. The interface BC was set between block five and block eight, and a crack with no contact between block two and block five was employed to show the advantages of this method in modeling and numerical simulations (Figure 10). The remaining settings were the same as in the previous example.

In this test, the wave field was more complex, which caused persistent disturbances by reflected and diffracted waves, but without spurious waves being generated by curve mesh conditions (Figure 11). We observed long-time disturbances and a reduced maximum of amplitude, which provided references for the seismic and postseismic analysis of the structure. This test also characterized the applicability and feasibility of SBP-SAT of the SMF of elastic wave equations to simulate structures with complex shapes.

6.4. Large-Scale Simulation for Complex Terrain

To demonstrate the full potential of the SBP-SAT methodology of the SMF in seismology, we next considered specific terrain in Lushan, China. Elevation data were derived from a profile of Baoxing seismic station, as shown in Figure 12, for an area 798 km in length. Since multiblock treatment involving complex terrain was the main focus of the numerical simulation, we employed homogeneous and isotropic parameters where  = 6400 m/s and  = 3695 m/s and blocks with grids discretized in each.

A method was developed combining PML [48] with the SBP-SAT methodology of the SMF of elastic wave equations, and a general curvilinear geometry was presented by this test. We focused on elastic wave propagation in the medium and the stability of the simulation. The initial data were set similar to that in (49).

The wave field at 5 s (Figure 13) showed that hypocentral P-wave propagated to the free surface, while Figure 14 showed the wave field at 25 s, once elastic waves had been reflected and scattered by the irregular surface, causing the wave field to be extremely complex. The wave field at 50 s and 70 s is shown in Figures 15 and 16, respectively. As can be seen, the numerical simulation led to a converging result. Since the PML methodology was used for the curve mesh, the absorption effect at ABCs was additionally strengthened and did not undermine the robust nature of the system.

Figures 15 and 16 show the wave field at 50 s and 70 s. The terrain was generally mountainous on the left and relatively flat on the right so that the wave field was more complex on the left due to the scattered wave. Figures 17 and 18 show the wave field at 100 s and 175 s. A long term and large-scale simulation with a stable wave field highlights the superior nature of using a multiblock approach in order to divide the large model into appropriate subblocks, which can improve simulation efficiency. Meanwhile, the continuous interface BCs had little impact on the wave field.

7. Conclusions and Future Work

The primary motivation of conducting this work was to derive a high-order accurate SBP-SAT approximation of the SMF of elastic wave equations, which can realize elastic wave propagation in a medium with complex cracks and boundaries. Based on the SMF, the stability of the algorithm can be proven more concisely through the energy method, which infers that the formula calculation and implementation of the program compilation for discretizing of the elastic wave equation will be more convenient, in addition to being able to expand the dimension of elastic waves. The SBP-SAT discretization of the SMF of elastic wave equations appears to be robust, despite the fact that an area is divided into multiple blocks, and the number of BCs is increased. The method developed in the present study was successfully applied to several simulations involving cracks and multiple blocks. It was also successfully applied to actual elevation data in Lushan, China, indicating that our method has the potential for simulating actual earthquakes.

In follow-up work, we aim to build on the current model with more accurate nonreflecting BCs, without undermining the model’s stability. Meanwhile, we also hope to include other SBP operators, such as upwind or staggered and upwind operators. Numerical simulation indicates that the method has wide application prospects. It can be challenging to model, discretize, and simulate actual elevation and geological structure, and strict requirements are necessary to ensure the stability of such a method.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 11872156), National Key Research and Development Program of China (Grant no. 2017YFC1500801), and the program for Innovative Research Team in China Earthquake Administration.