Joint Modeling of Wave Phenomena by Applying the Grid-Characteristic Method and the Discontinuous Galerkin Method

The aim of this work is to develop a hybrid computational method that combines the grid-characteristic method on regular structured grids with the discontinuous Galerkin method on unstructured tetrahedral grids. The proposed method makes it possible to describe integration domains with complex-shaped boundaries and contact boundaries and to calculate seismic fields taking into account the topography of the Earth’s surface, while saving computational resources. This modification of the method in the three-dimensional case has been proposed in this paper for the first time. Examples of using the developed method for calculating elastic wave phenomena arising during seismic prospecting are given. By way of testing, a comparison is made with results produced by the grid-characteristic method on curvilinear structured grids. The proposed hybrid method can be used not only for seismic prospecting, but also for calculating wave phenomena in other objects of complex shape.

These days hydrocarbon exploration requires the application of numerical methods for processing seismic survey data [1,2], including the development of high-accuracy methods for direct numerical simulation of wave phenomena in rocks.
In this paper, the grid-characteristic method on structured grids with constant step sizes in coordinates [10] is combined with the discontinuous Galerkin method on tetrahedral grids. Wherever possible, we use structured grids and the grid-characteristic method, while tetrahedral grids and the discontinuous Galerkin method are applied near boundaries and interfaces of curvilinear shapes.
Seismic wave propagation in rocks is simulated by solving the following initial-boundary value problem for the elastic wave equation: Here, is the velocity (derivative of the displacement), is a symmetric Cauchy stress tensor, r is the radius vector, is time, is the gradient vector, is the density, and are the respective speeds of P-and S-waves, is the unit tensor of the second rank, and denotes the tensor product of vectors: .
The free surface condition is set on the Earth's surface, and nonreflecting boundary conditions are specified on the other outer boundaries of the integration domain. As a source of seismic waves, we use a point source in the form of a Ricker wavelet with frequency 20 Hz.
Let us discuss in more detail the computation of the solution of system (1), (2) in the subdomains where the grid-characteristic and discontinuous Galerkin methods are used. In the case of the grid-characteristic method, the first step is dimensional splitting. Then, for each of three resulting one-dimensional hyperbolic systems of equations, we pass to Riemann invariants by applying the expressions Here, n 0 is the unit vector along the considered coordinate axis; n 1 and n 2 form a right-handed basis; ; and = .
Next, the resulting transport equations with coefficients are solved using their characteristic properties, and the back transition to the unknown functions and is made according to the expressions In the case of the discontinuous Galerkin method, the unknown functions are expanded in terms of a chosen basis of orthogonal functions in the unit tetrahedron : As a result, system (1), (2) is written in semidiscrete form: , , S ds (4) Here, are the coefficients of the solution expansion in terms of orthogonal functions in the tetrahedron m, where the first index corresponds to a component of the vector of unknowns; is the Jacobian for the coordinates of the nodes of the considered tetrahedron; and are outgoing and incoming fluxes; is the area of the jth face of the tetrahedron ; and the matrices , , and are expressed in terms of the coefficients of the hyperbolic system (1), (2) and the coordinates of the considered tetrahedron. The integrals in Eq. (4) are analytically computed in advance, since the set of orthogonal basis functions is known, and the resulting system is solved in time by an iterative method, for example, by Runge-Kutta methods.
Below, we describe the technique used to compute the solution near the boundary between the subdomains in which the discontinuous Galerkin method on tetrahedral grids and the grid-characteristic method on regular grids are applied.
Assuming that the tetrahedra along this boundary with a regular grid have identical orientations, we obtain the following simplified expression for incoming fluxes in boundary tetrahedra: Here, is the rotation matrix of the jth face; is the matrix made up of the eigenvectors of the hyperbolic system (1), (2); ,  To compute the unknown functions and by applying the grid-characteristic method, we need to know their values along characteristics lying in the tetrahedral grid zone. These values along characteristics are recovered using formula (3) in the boundary tetrahedron containing the considered characteristic.
Now we describe the results of test computations. Figures 1 and 2 show the three-dimensional geological model used in the test computations. Specifically, Fig. 1 presents the interfaces between the rock layers, and Fig. 2 depicts the distribution of the P-wave speed. The corresponding distributions of the other elastic parameters can be found in Table 1. Figure 3 shows the regular grids and the tetrahedral grids adjacent to the Earth's surface and to the curved interfaces between of the rock layers. The mesh size in the OX and OY directions was set to 12.5 m, and the       5 show the spatial distributions of the velocity modulus (wave patterns [14]) produced by the combined method and the grid-characteristic method on curvilinear structured grids [9], respectively. Tables 2 and 3 present the results of convergence analysis of the combined method proposed in this paper and the grid-characteristic method on curvilinear structured grids. The corresponding memory usage and the computation times are given in Table 4. By the step in coordinate, we mean the step size of regular grids in the combined method and the step size along the OX and OY axes in the method on curvilinear structured grids (since it was fixed in the considered formulation of the problem). Testing was conducted for a Courant number of 0.4. It can be seen that, for a sufficiently large step size, the proposed combined method does not provide advantages in terms of computational costs, but, under further mesh refinement, the fraction of the tetrahedral zone near the curved boundaries is reduced, and the computational cost is reduced as well. Additionally, the combined method has a higher order of convergence in both norms.
The results of the numerical experiments suggest that the developed numerical method is suitable for   Table 2. Results of convergence analysis of the grid-characteristic method on curvilinear structured grids [9] Step in coordinate, m high-accuracy solution of three-dimensional problems of seismic exploration of hydrocarbons and other deposits. According to the developed hybrid method, boundaries of complex geometry are described by constructing tetrahedral grids near them and applying the discontinuous Galerkin method, while the gridcharacteristic method on structured regular grids is used in the other subdomain of integration, which provides saving of computational costs. Note that the hybrid method can also be used to simulate wave phenomena in other objects of complex geometry.

FUNDING
This work was supported by the Russian Science Foundation, project no. 20-71-10028.

CONFLICT OF INTEREST
The authors declare that they have no conflicts of interest.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.