Optimized Mass and Stiffness Matrices for Accurate and Efficient 3D Acoustic Wave Modeling

Numerical modeling of the 3D acoustic wave equation regularly faces a trade‐off between accuracy and efficiency. Usually, an improved computational efficiency leads to a compromised numerical accuracy, in which the numerical dispersion error increases. Such a trade‐off is commonly reported for the finite difference methods, and the consistent, lumped, eclectic finite element matrices, and the high‐order spectral element methods. Therefore, we present a methodology that simultaneously optimizes both the accuracy and the computational efficiency of the finite element method for modeling the 3D acoustic wave equation in the frequency domain. The numerical dispersion objective function is derived and minimized for a tri‐linear interpolation function of a cube element in terms of the mass and the stiffness matrices. With four grid points per wavelength sampling, for less than 0.5% error in velocities, the optimized modeled solutions outperform the results from extant approaches. Numerical experiments with the homogeneous, 3D unbounded 2‐layer and the 3D Society of Exploration Geophysicists salt velocity models reveal that the optimized accuracy is persistent across all source‐receiver offsets, a performance observably absent from other finite element matrices. The optimized matrices are not only useful for numerical modeling, but are also readily available for accurate computation of the sensitivity matrices of the density and the bulk modulus parameters for the multi‐parameterization inversion.

Ideally, in principle we should seek a scheme that utilizes low order interpolation function with not-too-fine grids that can still efficiently provide good accuracy (i.e., least dispersion error). Therefore, based on the foregoing, we propose the optimization of the stiffness and mass matrices of the tri-linear interpolation function of the cube element. Our choice of the finite element method over the finite difference method is obvious. Like the finite difference method, the finite element is a numerical method for approximately solving the partial differential equations. Unlike the finite difference method, the benefits of the finite element method with the optimized matrices become all too clear in the multi-parameterization inversion, that is, where the simultaneous computation of the sensitivity matrices for the density and the bulk modulus parameters can be made. Our methodology involves the minimization of an objective function of the normalized phase velocity for the purpose of conducting a dispersion analysis. Illustrative velocity model examples will be subsequently cited to quantify the accuracy of the optimized matrices relative to the lumped, consistent and eclectic matrices presented in Marfurt (1984). The velocity models that will be used are the 3D homogeneous, the 3D unbounded two-layer (to validate the accuracy of the reflection energy), and the 3D SEG (Society of Exploration Geophysicists) salt velocity models.

Finite Element Stiffness and Matrices of Tri-Linear Function
The source-free 3D acoustic wave equation in the frequency-domain is written as: where the 3D pressure wave field, p = p(x,y,z,ω), can be modeled for a velocity model V 0 at the angular frequency of ω. By variational calculus, the solution p in Equation 1 in the computational domain D, is such that will also minimize the function (Norrie & de Vries, 1978): Domain D can be divided into n finite elements in order to sum the contributions of each element to the total χ. Thus, each element functional χ e is from the element trial function , whose differentiation with respect to the element nodal vector p e gives: where the element stiffness matrix K e and element mass matrix M e are given by where N i ,N j are the tri-linear shape functions and i,j identify node 1, 2, … of the elements.
It is instructive to note that the element stiffness matrix can be expanded as: Using the tri-linear shape function on a cube element of a side length h, the element stiffness and mass matrices are: where a = 4; b = 2; c = 1; d = 8; e = 4; f = 2; g = 1.

Derivation and Minimization of Objective Function
The inner node (i,j,k) of a cube element in Figure 1 is surrounded by 26 nodes. Assembling the stiffness and mass matrices, from Equation 7 to Equation 10, at the inner node, leads to the values in Table 1 for the inner node and its 26 neighbors. The 3D plane wave equation is: where ι = √ −1 ,k x = Γ sin θ cos ψ,k y = Γ sin θ sin ψ,k z = Γ cos θ are the wavenumbers in the x,y and z directions respectively; Γ = ω/V o ; θ is the propagation angle from the z-axis and ψ is the azimuthal angle.
Assembly of the element matrix equations (see values in Table 1) yield a system of equation of the form: Now from Equation 12, we can express the wavenumber Γ = ω/V 0 as: where V ph is the phase velocity.
From Equation 23, we have Or, where the number of grid points per wavelength (λ) is G = λ/h = 2π/Γh.
The objective function to be minimized is thus given as: We used the numerical method presented in Jo et al. (1996) to minimize the objective function in Equation 26. Table 2 shows the optimized matrices inverted from Equation 26, as well as those for the consistent, lumped and eclectic. The normalized phase velocity plots in Figure 2 (using the display software from She, 2019) demonstrate the dispersive effects of the different matrices. The numerical anisotropy values are 9%, 15%, 6% and 0.15% for the consistent, lumped, eclectic and the optimized matrices respectively. While the optimized matrices ( Figure 2d) produce the least dispersive effect on the modeled acoustic wave field, the dispersion associated with the lumped matrices is the largest, increasing as the grid size becomes bigger (Figure 2b). Further analysis indicates that the number of grid points per wavelength to produce 1% error in velocities, for the eclectic, consistent and the lumped models are respectively 6.21, 12.3 and 14.1. Comparatively, 4 grid points per wavelength unit sampling are required to produce 0.45% error in velocities, factors that will simultaneously improve both the computational efficiency and the numerical accuracy of the modeled acoustic wave field. Next, we will perform some numerical experiments with the different matrices to qualitatively evaluate their dispersive effects on the generated seismograms.

Homogeneous Velocity Model
With a point source located at (x s ,y s ,z s ), the analytic 3D Green's function, in the frequency-domain, in a homogeneous velocity medium can be computed with (Arfken & Weber, 2005):   Results from using the different matrices are compared with the analytical solution in Figure 3. Trace 70, from all the matrices, is overlaid on the analytical trace for comparison in Figure 3a. For the lumped matrices, the computational efficiency is favored as numerical accuracy becomes highly compromised; and evidently, the lumped trace demonstrates the full effects of dispersion error on the modeled seismogram. These effects are more pronounced after the first arrival energy. Consistent and eclectic results are comparatively the same, with subtle noticeable dispersion error before the first arrival energy. Moreover, their normalized energy amplitudes are observably smaller than the true solution. Unlike the normalized traces from other matrices, the normalized optimized trace rids the energy of dispersion error, with accompanying higher resolution than others. Figure 3b presents the performance of the optimized matrices at both near and far offsets. The reductions in the dispersion error and increased resolution are maintained even as source-receiver offset increases. Moreover, for the lumped, consistent and eclectic matrices, it is noteworthy that the dispersion error and amplitude inaccuracies increase with offset.

3D Unbounded 2-Layer Velocity Model
In a 3D unbounded 2-layer velocity model, a horizontal surface separates two layers which are acoustically different. Layer velocities and densities are respectively given as c n and ρ n , where n = 1,2. The closed-form solution for such a problem is given in Officer (1958). Therefore, this velocity model will facilitate the validation of the accuracy of the optimized matrices on reflected energy, which is apparently absent in the homogeneous velocity model example.
We use the layer velocities of (c 1 ,c 2 ) = (1500,2500)m/s and keep their densities fixed at 1,000 kg/m 3 . The first derivative Gaussian source pulse, with maximum frequency of 15 Hz, is positioned at (x s , y s , z s ) = (1,250,1,250,500) m, and the receiver is buried at the depth of 500 m. The number of elements in the x, y, and z directions are respectively 300, 100 and 100, with a constant grid size of 25 m. Similar to the homogeneous velocity model case, the frequency domain solutions are produced for 61 frequencies as the PML boundary condition is implemented.
With the aid of Figure 4, we compare the optimized seismogram with the analytical solution. Figure 4a compares the accuracy of the analytical solution, for trace 70, with the different mass and stiffness matrices. The direct wave and reflected energy are captured in trace 70. Again, the lumped trace significantly deviates from the true trace, while the amplitudes of the direct and reflected energies of both the consistent and eclectic traces visibly differ from the true trace. However, the optimized trace matches appreciably well with the analytical trace. In the near offset, the amplitude of the direct wave energy dominates the seismogram, and the optimized traces agree with the true traces as seen in Figure 4b. Parts of the direct wave energy, as well as parts of the reflected energy are present in the mid-offset, where there is no lack of fit between the optimized traces and analytical traces; see Figure 4c. In the far offset, the direct wave and the reflected wave energy have constructively interfered, and yet the accuracy explanation carries through the same reasoning as for the other offsets. We note that although the optimized matrices maintain accuracy of the results across all offsets, the same conclusion cannot be said of the other matrices, whose accuracies deteriorate with offsets.

SEG Salt Velocity Model
The SEG salt velocity model in Figure 5a is 10 by 10 by 4 km 3 in extent, and it is discretized by cube elements, with side length of 60 m each. The SEG salt velocity model is a highly faulted, and laterally extensive high velocity salt body; the number of elements in x, y and z directions are 169, 169 and 67 respectively. The source wavelet of maximum frequency of 7 Hz, is planted at (x s , y s , z s ) = (6780, 6780, 100) m, while the receivers are placed at the depth of 100 m. Imposing the top Dirichlet boundary condition, the inverse Fourier transform of 36 wave fields in the frequency domain leads to the computation of an optimized 3D shot gather. The shot gather is presented in three different views namely: the top at 0 s in Figure 5b, the top at 1 s in Figure 5c, and the top at 4 s, inline at 6.6 km and crossline at 6.78 km. Indeed, with the optimized matrices, high resolution of the reflected,

Discussions
The use of the optimized 3D tri-linear stiffness and mass matrices, from the minimization of dispersion objective function, has proved to be both accurate and efficient. This is not so for the lumped, consistent, and the eclectic matrices that trade-off accuracy for efficiency and vice versa. With the 4 grid points per wavelength unit sampling required to produce 0.45% error in velocities, the optimized matrices do not only outperform the published 4.5 grid points per wavelength unit sampling of the high-order spectral element method (SEM), but also return higher resolution images than the SEM. In comparison with the finite difference method, the optimized matrices offer accurate multi-parameterization inversion, that conveniently does simultaneous computation of the sensitivity matrices for the density and the bulk modulus.
On the last note, the computation time (t) for the lumped, consistent, eclectic and optimized methods can be defined as a function of the bandwidth (B), which is a function of number of points per wavelength: where B = 2q, q is the number of grid points per wavelength and s is the spatial dimension, which 3 for our 3D modeling experiment. Substituting the values of the number of points per wavelength for both eclectic, 6.21, and the optimized, 4, methods into Equation 28, the optimized method is 1.56 smaller in bandwidth than the eclectic method. Consequently, the optimized method is at least 3.7 faster than the eclectic method. Obviously, the computational efficiencies for the lumped and consistent are poorer than eclectic since they have bigger bandwidths.

Conclusions
The mass and stiffness matrices of the tri-linear interpolation function of a cube element have been optimized by minimizing the normalized phase velocity objective function. Consequently, the optimized matrices efficiently model accurate 3D acoustic wave field. Accuracies of the transmitted and reflected energies in the homogeneous and the 2-layer velocity models for different matrices are compared with the analytical solutions. The optimized results outrank solutions produced with the consistent, lumped and eclectic matrices in accuracy and resolution. Without fail, the optimized accuracy is remarkably maintained across all offsets. The 4 grid points per wavelength unit sampling required by the optimized matrices also outclass the published results of the spectral element method. It is therefore envisaged that the optimized matrices will be useful in both acoustic wave modeling and in inversion, particularly the multi-parameterization inversion of the density and the bulk modulus.

Data Availability Statement
Our numerical models do not utilize (input) data, but the figures are already available in electronic form. The 3d Figure