Computer Simulation of the Seismic Wave Propagation in Poroelastic Medium

: This article presents an algorithm for the numerical solution of an initial ‐ boundary value prob ‐ lem for a symmetric t ‐ hyperbolic system of partial differential equations. This problem is based on con ‐ tinual filtration model, which describes the propagation of seismic waves in a poroelastic medium satu ‐ rated with a fluid characterized by such physical parameters as the propagation velocities of longitudinal P ‐ (fast and slow) and transverse S ‐ waves, the density of the medium materials, and porosity. The system of linearized equations of saturated porous media is formulated in terms of physical variables of the ve ‐ locity–stress tensor of the porous matrix and the velocity–pressure of the saturating fluid in the absence of energy dissipation. The solution is implemented numerically using an explicit finite difference upwind scheme built on a staggered grid to avoid the appearance of oscillations in the solution functions. The program code implementing parallel computing is developed in the high ‐ performance Julia program ‐ ming language. The possibility of using the approach is demonstrated by the example of solving the problem of propagation of seismic waves from a source located in the formation. Computational experi ‐ ments based on real data from oil reservoirs have been implemented, and dynamic visualization of so ‐ lutions consistent with the first waves arrival times has been obtained.


Introduction
Mathematical modeling plays an essential role in the development of new geophysical technologies, which makes it possible to provide synthetic wave fields for various arrangements of sources and receivers of acoustic vibrations in models of media with a given complexity of internal structure. Computer modeling is especially important in view of the complexity of conducting an experimental study of internal non-stationary physical processes in fluid-saturated elastically deformable porous structures [1,2]. The development of computer simulation methods for realistic models of filtration mechanics in geophysical applications stimulates progress in many other fields of research, such as medicine, biotechnology, materials science, chemistry, micro-and nanofluidics, and geosciences [3][4][5][6][7]. Hyperbolic systems of partial differential equations are found in modeling separation and reaction engineering problems.
The main approach to the study of wave processes in liquid-saturated porous media is based on the Frenkel-Biot theory [8,9]. The model, which describes the processes of deformation of an elastic porous medium and the flow of fluid in it, is macroscopic, i.e., assumes that the space is filled with an enclosing poroelastic two-phase medium, and the phases corresponding to the porous solid and the liquid contained in the pores are present at each point in space. The Frenkel-Biot type models consist of two groups of equations: (I) equations of the theory of elasticity taking into account the internal forces associated with the influence of fluid pressure in pores on the stress-strain state of the medium, and (II) equations of fluid filtration in pores taking into account changes in the volume of the pore space filled with fluid due to deformation of the medium.
A continuum filtration model that satisfies the laws of thermodynamics was proposed in [10,11]. In this theory, the properties of a saturated porous medium are described using such physical parameters as the propagation velocities of two longitudinal waves and one transverse acoustic wave, the partial densities of the porous matrix and the saturating liquid, and the permeability coefficient. The continuum filtration approach makes it possible to describe the poroelastic medium with a smaller number of model parameters than in the Frenkel-Biot theory of poroelasticity and with a better agreement with experimental data [12].
The main tool for studying the applied problems of the theory of poroelasticity is numerical modeling, which allows conduction of numerical experiments with fully controlled input data and thereby reliably verifying hypotheses and formulating criteria for the manifestation of fluid mobility in seismic data. In [13][14][15][16][17][18][19][20][21][22][23][24][25], the propagation of seismic waves in porous media in the presence of energy dissipation was studied. Finite difference methods for solving poroelasticity problems have been formulated in several ways. In [13,14], the central difference method was used in terms of displacements. A semi-analytical method for solving a system of poroelasticity equations in terms of displacements is proposed in [15,16]. In [17], the problem was investigated by means of a predictor-corrector scheme for a system of velocity-stress equations. In [20], a spectral method was used, where the time derivatives are approximated through Laguerre functions (polynomials). The numerical solution of linear two-dimensional dynamic problems for porous media by the finite difference method together with the discrete PML model was used in [19] for the absorbing boundaries. The PML is introduced by Berenger (1994) for electromagnetism. Many authors have applied the PML to various wavefields simulations: Zeng et al. (2001) analyzed the PML in numerical modeling of the poroelastic media in time domain [14]. In [2], the PML is applied to the poroelastic waves in frequency domain and is used for the absorbing boundaries. The simplified wave equation with the PML attenuation term in time domain is used in [3]. Computational algorithms based on the finite element method have been successfully used to simulate the dynamics of hydraulic fracturing crack development in a three-dimensional formulation [26]. The development of multiprocessor systems stimulated the development of computational methods for solving seismic problems [27,28]. The grid-characteristic numerical method, which takes into account the internal mathematical structure of the hyperbolic problem-the propagation of discontinuities along the characteristics-has been successfully applied to the numerical solution of the direct problem of wave propagation in a heterophase medium [29].
In this paper, a numerical solution of the dynamic poroelasticity problem is obtained within the framework of a continuous filtration model. The features of the discretization of linearized poroelasticity equations of physical variables are considered. A feature of the model of a homogeneous, isotropic liquid-saturated porous medium is the specification of its properties by experimentally measurable velocities of longitudinal fast and slow P and transverse S acoustic waves in a two-phase medium. The difference algorithm is based on the finite difference method using an upwind staggered difference scheme to avoid the appearance of oscillations in solution functions [30]. The implementation of the scheme is sufficiently universal for calculations on a rectangular grid with the setting of various boundary conditions for pressure and stresses. The program code for the implementation of the algorithm was developed using the high-performance Julia programming language, since the explicit computational scheme allows efficient parallelization of cyclic sections of the algorithm. The stability of the scheme is ensured by selecting the grid steps in time and in spatial variables to satisfy the Friedrichs-Courant-Lewy criterion. The possibilities of the approach used are demonstrated by two examples of numerical experiments using realistic sets of parameter values and types of functions describing a seismic source. A visualization of the solution describing the change in the velocity of propagation of seismic waves in a poroelastic formation is obtained.

Statement of the Problem
Let us turn to the mathematical formulation of the continuum model for a two-dimensional dynamic problem. It is expressed as the initial-boundary problem for a symmetric t-hyperbolic system. The main equations are based on the laws of momentum conservation, Hooke's laws, and are consistent with the thermodynamics conditions. A halfplane filled with a liquid-saturated porous medium is considered, with parameters characterizing each of the components of such a medium.
It is assumed that the liquid completely fills the pores. The distribution of phases in space is described by such a macroscopic parameter as porosity (the volume concentration of voids filled with fluid in the medium).
The propagation of seismic waves in these media in the absence of energy loss is described by the following initial-boundary value problem in terms of the velocity components of the saturating fluid, the velocity of the solid matrix, the fluid pressure, and the stress tensor [19][20][21][22]24,25]. The values related to the solid phase of the elastic medium are denoted by the subscript "s", and to the mobile phase of the saturating liquid-by the subscript "l". Thus, partial density of an elastic porous medium is denoted as s  and a partial density of a liquid is denoted as l  . Partial densities are determined through the physical densities of the porous medium f s  and liquid f l  , and porosity-, and together they form a density of the porous medium 0  , as follows: The problem considered in the domain stated as the following system of equation: The equation of motion for an elastic medium: where for the two-dimensional case 1, 2 j  , is expressed in the form of conservation laws in the linearized case, where is the velocity vector of a liquid; Hooke's law for a solid matrix (elastic medium): where jk  is the Kronecker symbol; Hooke's law for a saturating liquid is expressed as follows: Initial conditions: Boundary conditions on the free surface in the plane: In the two-dimensional case, the Equations of system (1)-(4) consist of eight partial differential equations of the following form: Boundary condition (6) only for a porous layer saturated with a liquid with a density 0 l   takes the following form with respect to the functions 22  and p : We write the system (7)- (14) in a matrix for 12 22 , , , , , , , To find a numerical solution of the dynamic problem for the poroelasticity Equations (7)-(15), a discrete analogue is constructed. An algorithm for computing a solution based on an explicit scheme is built on a uniform staggered grid using the finite difference method.
Since in problem (7)-(15) the boundary conditions are set only for three functions (15) on the left boundary at 2 0 x  , for the remaining functions on the right boundary of a sufficiently large area, the conditions for the decay of functions at infinity are assumed.
To obtain the numerical solution, we use an upwind difference scheme, which approximates the derivatives, starting from the right boundary (a rather large value is taken instead of infinity) along the 2 x -axis moving to the left to 2 0 x  . To ensure the convergence of a numerical solution based on an explicit difference scheme, it is necessary to choose steps in time and space variables so that they satisfy the Courant-Friedrichs-Lewy (CFL) conditions, i.e., the following inequalities hold: x axis, and 2 h is the step along 1 x axis.
The flowchart of the numerical solution is shown on Figure 1:

Discrete Approximation of the Problem
The use of the finite difference method has shown its effectiveness for the numerical solution of the problem stated on the constructed uniform difference grid in the region is specified with a staggered arrangement of nodes. Let us introduce the following notation for such a uniform difference grid:   22  the stress tensor and pressure are in Due to the presence of numerous indices, we introduce the following designations for the desired grid functions with the same letters as in the differential problem (7) , where coefficients i  , are calculated by the following formulas: The calculations were implemented using parallel computing on a computer with the technical characteristics shown in Table 1:

Description
Technical Characteristics Processor 16-core AMD Ryzen 9 3950X Clock speed/Frequency 3.5 GHz (Matisse) RAM 64 GB The use of an explicit computational scheme, which contain the cyclic sections, enables the software to create parallelized algorithm and allows the efficient use of the computational resources.
The dynamic visualization of the computational experiments results presented in this paper was performed using ParaView software package.

First Computational Experiment
The real input data for physical parameters in the first computational experiment of the problem are taken from [19] and are shown in Table 2: Table 2. Input physical parameters and grid characteristics in the first computational experiment. Step along axis 1 Step along axis 2 Step in time The source function by the following formulas: where the momentum function is described by the formula: 14).
The graph of the function ) (t f is shown in Figure 4.  On the dynamic graphs on Figures 5 and 6 we can observe that the velocity in a solid phase change faster in time that in a liquid phase, which is in agreement with the observed experimental data. The seismograms of components 1 2 ( , ) u x t and 2 2 ( , ) u x t along the traces going through the source point are shown in Figures 7 and 8, respectively.

The Second Computational Experiment
In the second computational experiment, the different set of values of physical parameters shown in Table 3 and the form of the source function describing the seismic pulse of the problem are taken from [20]: Table 3. Input physical parameters and grid characteristics in the second computational experiment. Step in time

##
The delta function   0   x x and its partial derivatives are approximated in a similar way as in the first computational experiment using the "hat" function. The graph of the impulse function ) (t f is shown in Figure 9.   The program code implementing parallel computing is developed in the high-performance Julia programming language. The possibility of using the approach is demonstrated by the example of solving the problem of propagation of seismic waves from a source located in the formation. Computational experiments based on realistic data from oil reservoirs have been implemented and dynamic visualization of solutions consistent with the first waves entry times has been obtained.

Discussion
Thus, an original algorithm for the numerical solution of the initial-boundary value problem for the poroelasticity system is presented. This problem models the propagation of seismic waves in a poroelastic medium saturated with a fluid, characterized by such physical parameters as the propagation velocities of longitudinal and transverse waves, the physical density of the medium materials, and porosity. This approach has important industrial applications. Modeling of seismic waves propagation is instrumental for seis-mic surveys as the most reliable geophysical technique used to identify oil and gas prospects in the geologic structures. This paper introduces a new approach to the computer simulation of wave propagation processes in complex multiphase media. The proposed algorithm allows efficient use of parallel computations on computers with multi-core processors. Thanks to improvements in high-performance modern computing systems, the introduction of a new approach to computer simulation technology may play a critical role in the field of data processing and interpretation and in other applications. The use of a staggered grid method enabled us to avoid oscillations of the eight solution functions in a two-dimensional case of the problem. The solution is implemented using an explicit finite difference upwind scheme built on a staggered grid using parallel computing. The continuum filtration approach makes it possible to describe the poroelastic medium with a smaller number of model parameters than in the Frenkel-Biot theory of poroelasticity and with a better agreement with experimental data.
The numerical modeling allows the conduction of numerical experiments with fully controlled input data and thereby reliably verifying hypotheses and formulating criteria for the manifestation of fluid mobility in seismic data. The program code implementing parallel computing is developed in the high-performance Julia programming language. The stability of the proposed scheme is ensured by the satisfaction of the Friedrichs-Courant-Lewy criterion. Computational experiments based on real data from oil reservoirs have been implemented, and dynamic visualization of solutions consistent with the first waves arrival times has been obtained.
This work is mainly aimed at demonstrating the numerical method for solving the stated poroelasticity problem, while the results of the theoretical study of the resolvability conditions and the expression of the analytical solution in explicit form will be the subject of a separate paper.
This investigation presents an algorithm for the numerical solution of an initialboundary value problem for a system of symmetric partial differential equations of t-hyperbolic type. This will be crucial for research in various areas, such as medicine, chemical engineering, micro and nanofluidics.

Patents
The authors are planning to create an easily usable interface allowing implementation of solutions for various input data of physical parameters values and a source function form based on the algorithm described in the article. Once this work is completed, a patent will be obtained resulting from the work reported in this manuscript.