Simulation of 2D Acoustic Wave Propagation with Absorbing Boundary Condition

In this work we design an algorithm to simulate the 2D acoustic wave (P wave) propagation in a medium with different regional velocities. We use the Generalized Conjugate Residual (GCR) algorithm to solve the trapezoidal time integration, where for a given time step of 0.01 seconds, the error is approximately 8% with respect to the reference. We further speed up the simulation by a factor of 900 by implementing a model reduction scheme using Proper Orthogonal Decomposition (POD). Lastly, to suppress the artificial reflections formed when the wave bounces off the fixed boundaries, we created two virtual absorbing layers by tapering the solution of nodes close to the corresponding boundaries after each time iteration. However, we note that the successful implementations of the model reduction and absorbing boundary condition are not compatible with each other and provide explanations for why and how to combine them in future work.


Introduction
Earthquakes generate seismic energy waves that propagate the earth's interior and get recorded by geophones when they reach the surface, as shown in Fig. 1. The seismic recordings are then analysed to extract geological information about the Earth. The simulation of seismic wave propagation is widely used in the seismology community for various purposes. Sometimes, the synthetic data generated by simulation is treated as the theoretical replacement of real data to be tested with novel imaging processing techniques. For instance, techniques like the curvelet transform and slant stacking are used on the seismic data to analyse the Earth's mantle [1,2]; however, the real data are usually contaminated with ambient noise and false instrumental response. Thus, synthetic data is first used before any work on the real data. Since different earthquakes generate different seismic data, seismologists wish to generate simulations for different scenarios as quickly and as accurately as possible. To cater to these demands, we focus on speeding up the simulation and mimicking the real-world scenarios by adopting model reduction and absorbing boundary conditions. However, we simplify the Earth's geometry to a rectangular 2D grid due to the limited scope of this work. [3] Figure 1. Illustration of P-wave propagation and the velocity model.

Problem formulation
Our work starts from the 1D wave equation: where ≡ . The function is synthesized into a linear system, where is made of the finite difference discretization matrix and identity matrix : We setup the grid as a × 2 network, where is the number of nodes.
• The state of our system is a vector of length 2 and each node includes displacement and particle velocity : • The parameters of our components includes: where is the P-wave velocity, ∆ is the spatial increment, is the model distance, is the location of excitation, and ∆ is the time increment.
• The excitation of the system is: • Thus, the final set of equations are written explicitly as: • The output is the displacement of surface nodes as a function of both distance × ∆ and time , which is known as the synthetic seismic recording at the surface.

Model parameters
In this study, ∆ = 200 , ∆ = 0.01 . = 6600 / for the top half block and = 9900 / for the bottom half block. The total horizontal and vertical distances are 18000m and 9000m respectively, which yield a total number of 4186 simulating points. The acoustic wave velocity is the physical property assigned to each node as constant. The scope of this study is to simulate the first 5 seconds of the wave propagation.

ODE integration method
We compared the Forward Euler (FE) and Trapezoidal approximation methods, and found that Trapezoidal is preferred. During the FE progression, the system diverges easily because of the accumulated errors, which resemble introducing fake energy into the system. Thus, the largest feasible ∆ for FE is 10 −6 , which is computationally expensive. Trapezoidal on the other hand, increases ∆ to 10 −2 . The structure of algorithm is shown in Fig. 3. We use GCR iteration inside Newton method to solve for each time interval. Note that, the system is completely linear and does not require 6 Newton method, but it is retained for the future implementation of nonlinear terms, such as energy loss. In our case, GCR takes about 3.3 min to finish, compared to around 23 min for LU decomposition.

The technical challenge
We completed two technical challenges: model order reduction using the Proper Orthogonal Decomposition (POD) method and absorbing boundary conditions.

POD model order reduction
In seismology, only nodes at the surface matter because that is where geophones are located. Therefore, model reduction is desired to speed up the simulation by preserving only the input-output relationship between the excitation and surface nodes. With POD, singular value decomposition (SVD) is applied to matrix , constructed with a sufficient number of 's generated by the original model (501 in this case). Projection matrix is constructed by taking the singular vectors corresponding to largest singular values, in our case = 140. The model can then be reduced with ̂, ̂ and ̂ and can quickly and effectively handle different excitations as shown in Figure 4.      Fig. 7 shows a comparison between the original model and the reduced model. Both look very similar, and specifically the −2 error between the two is only 0.0124%. However, the reduced model only takes 0.22s to run, while the original one takes about 25 min, a speed up factor of 900.

Technical discussion
In Fig. 8, the reciprocal of the slope where the hyperbola converges is 6000 m/s for the top and 9000 m/s for the bottom, which is close to the true P-wave velocities assigned to the top and bottom blocks ( = 6600 / , = 9900 / ). The velocity approximations together with the half-depth reflection suggest that our model is working properly, and the results are valid. Although our reduced model makes the simulation 900 times faster than the original model, one can easily speed it up more. In this study, we pick 140 largest singular values, and achieve a low error of 0.0124%, which may not be necessary in other applications. Thus, choosing fewer singular values by making smaller is feasible to further speed up the simulation. Since we only take 10% of nodes from both sides to be the absorbing layer, the result still experiences faint side reflections, shown in 9. If one increases the size of the layer, the performance of suppression will be better, but it will become more computationally expensive. Implementing the absorbing boundary condition using tapering window has its limit. The tapering process after each time iteration is manipulating matrix in GCR, but the model reduction requires to be unchanged throughout the simulation. Thus, our implementation of absorbing boundary condition is not compatible with the POD model reduction. The address it, more advanced implementation of absorbing boundary condition is needed. For instance, Neumann boundary condition is working directly on matrix , but it requires one to consider first order spatial derivative, which adds another layer of constrains to the model [4][5][6][7].