{'fraction_non_alphanumeric': 0.03361062799907801, 'fraction_numerical': 0.018670242859836976, 'mean_word_length': 4.373719175768494, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Two-phase hyperbolic model for porous media saturated with a viscous fluid and its application to wavefields simulation

We derive and study a new hyperbolic two-phase model of a porous deformable medium saturated by a viscous fluid. The governing equations of the model are derived in the framework of Symmetric Hyperbolic Thermodynamically Compatible (SHTC) systems and by generalizing the unified hyperbolic model of continuum fluid and solid mechanics. Similarly to the unified model, the presented model takes into account the viscosity of the saturating fluid through a hyperbolic reformulation. The model accounts for such dissipative mechanisms as interfacial friction and viscous dissipation of the saturated fluid. Using the presented nonlinear finite-strain SHTC model, the governing equations for the propagation of small-amplitude waves in a porous medium saturated with a viscous fluid are derived. As in the conventional Biot theory of porous media, three types of waves can be found: fast and slow compression waves and shear waves. It turns out that the shear wave attenuates rapidly due to the viscosity of the saturating fluid, and this wave is difficult to see in typical test cases. However, some test cases are presented in which shear waves can be observed in the vicinity of interfaces between regions with different porosity.


Introduction
Modeling processes in saturated elastic porous media is of interest to many areas of industry and our daily life. The research in this area began in the middle of the last century in the works of Biot [3,2]. In the subsequent development of Biot's theory, special attention was paid to seismic applications modeling, such as small amplitude wave propagation and energy loss, see, for example [4,13] and references therein. The main feature of wave propagation in Biot's model, in comparison with the theory of elastic waves, is that there are two compressional waves in it -fast and slow. Herewith, the slow compression wave attenuates quickly and therefore it is difficult to detect it (it can be clearly seen only for the high frequencies). The attenuation mechanism of this wave is an interfacial friction, which is also an important factor in the loss of wave energy in a porous saturated medium.
In Biot's theory, there are no assumptions about the rheology of the saturating fluid and an interfacial friction is taken into account as a result of the average motion of the fluid relative to the solid skeleton [13]. Fluid viscosity is indirectly presented in the friction coefficient, because its definition is based on the Darcy law (limiting case of Biot's theory). To the best of our knowledge only in a few papers the study of influence of saturating liquid rheology on wavefields can be found. Of particular note is the paper [12] that treats heavy oil as a viscoplastic saturating fluid, which is not of interest for the present paper. There is an attempt to include the viscous stress directly into the model of the saturating fluid, which is assumed to be Newtonian [26,14,8] (see also the book [13]). The effect of viscosity on the wavefield leads to the appearance of the so-called slow S-wave. As it is noted in [26], this wave decays rapidly, therefore it is difficult to detect it even in numerical simulations. In [26] it is reported that S-wave can be found only in the vicinity of interfaces. Thus, taking viscous stresses in a saturating fluid into account provides an additional mechanism for energy loss in a porous medium and should be taken into account in the interpretation of seismic data.
In fact, a saturated porous medium can be considered as a two-phase mixture of a saturating fluid and a solid skeleton. This means that methods of the continuum mechanics theory of multiphase media can be used. In papers by Wilmanski [28,27], the correspondence between Biot's theory and the theory of a two-phase fluid-solid mixture was investigated. From these papers, it can be concluded that Biot's model in its governing equations can be interpreted as a continuous two-phase model under certain assumptions.
Recently, a new Symmetric Hyperbolic Thermodynamically Compatible (SHTC) model of a two-phase continuum has been proposed for an elastic porous medium saturated with a compressible fluid [23]. The SHTC system theory allows formulating a class of symmetric hyperbolic systems of partial differential equations that satisfy the laws of thermodynamics and have good mathematical properties [11,10,9,21,24,16]. This class includes the governing equations of many known models of continuum mechanics, and moreover SHTC systems theory can be a powerful tool for developing new models of complex processes. The finite-strain poroelastic model presented in [23] is based on the unified model of continuum [17,6,7] which is capable to describe a medium in elastic, elastoplastic and viscous fluid states using a single system of governing differential equations. It extends the unified model by considering a saturated porous medium as a two-phase solid-fluid mixture, and its governing equations combine the equations of the unified model and equations of the SHTC model for compressible two-phase fluid flows [20,19,22].
In [23], the rheology of the saturating fluid is not taken into account, and the only mechanism of energy loss in wave processes is interfacial friction. In the present paper we extend the model [23] by introducing the viscosity of the saturating fluid in the model. This extension is based on the use of a hyperbolic viscous fluid model, which is essentially part of the unified model [17,6] and represents a viscous fluid as a solid medium with very small shear stress relaxation time. Thus, the two-phase continuum model of porous medium discussed here is considered as a mixture of two viscoelastic constituente, one of which corresponds to the fluid. The closing relations (elastic moduli and mixture shear stress relaxation time) in the model are taken as functions of the volume and mass fractions of constituents and determined by some mixture rules. This allows one to obtain the governing equations for a pure solid and a pure fluid as limiting cases of the model system with vanishing volume fractions of the constituents. It is important that the presented model is hyperbolic and guaranties finite velocities for all types of propagating waves.
Limiting the full non-linear finite-strain model to the case of small amplitude wave propagation, we derive a linear poroelastic model that can be used to study wavefields in a porous medium saturated with a viscous fluid with inclusions of pure solid and pure fluid regions. This is very advantageous property because it allows one to simulate waves in regions with a complex geometry of the internal structure by the single partial differential equation (PDE) system. For example, in [18] the so-called diffuse interface approach is used to simulate wavefields in porous media with pure solid and pure fluid inclusions, which demonstrates the effectiveness of the model.
The results of the analysis of the propagation of one-dimensional waves show that there are three types of waves in the model: fast and slow compressional waves and shear waves (as in Biot's theory). The appropriate choice of closing relations gives a good qualitative agreement of the frequency-dependence of both compressional waves with those obtained by Biot's theory. The behavior of the shear wave strongly depends on the parameters that determine the effective viscosity of the medium: the velocity of the shear wave tends to zero at low frequencies, and the corresponding attenuation factor increases with increasing viscosity. What concerns the observation of the so-called slow S-wave, in the model there is only one shear wave, but in some circumstances its velocity can be small. This results in that, in numerical simulations, one can see some smeared waves at the vicinity of the interfaces between media with different porosities, which is similar to the observations given in [26]. Thus, the presented model is capable of simulating wave phenomena in a porous medium saturated with a viscous fluid. Yet, for a correct quantitative description of shear wavefields, reliable experimental data on their propagation are needed. The rest of the paper is organized as follows. In Section 2, the general two-phase solid-fluid SHTC model is formulated in the finite-strain settings and its properties are studied. Section 3 contains a description of a simplified isentropic single pressure model and its limiting versions for vanishing phase volume fractions. In Section 4, a derivation of the equations for small amplitude wave propagation in a porous medium saturated with a viscous fluid is presented. In Sections 5 and 6, we analyze characteristic velocities and dispersion relations (frequency-dependence of the sound speeds) of the model. Finally, in Section 7, a series of two-dimensional test problems is solved using a finite difference scheme on staggered grids. In particular, an observation of the slow S-wave is discussed.
2. Two-phase thermodynamically compatible model of solid-fluid mixture applicable for the description of saturated porous medium In this section, we consider the two-phase solid-fluid mixture model [23], also called master model here, that later will be used for the derivation of PDEs for small amplitude wave propagation is based on the model for deformed porous medium saturated with the compressible fluid. If solid and fluid constituents are characterized by their volume fractions α 1 and α 2 , (α 1 + α 2 = 1), then the master system reads as Here, ρ = α 1 ρ 1 + α 2 ρ 2 is the mixture density, ρ 1 , ρ 2 are the fluid and solid mass densities, c 1 = α 1 ρ 1 /ρ is the fluid mass fraction (c 2 = 1 − c 1 = α 2 ρ 2 /ρ is the solid phase mass fraction), v i = c 1 v i 1 + c 2 v i 2 is the mixture velocity, v i 1 , v i 2 are fluid and solid phase velocities, w i = v i 1 − v i 2 is the relative velocity of phase motion, s is the entropy of the mixture, A i k is the distortion matrix characterizing local elastic deformation of the mixture. Equation (1a) is the mixture momentum equation, (1b) is the evolution equation for the elastic distortion of the mixture A i k , (1c) is the mixture mass conservation law, (1d) is the fluid mass conservation law, (1e) is the balance equation for the relative velocity, (1f) is the balance equation for the liquid volume fraction, and (1g) is the mixture entropy balance law. The source term in the latter equation provides the growth of the mixture entropy due to the dissipative processes as the second law of thermodynamics requires.
The most important closing relation for system (1) is the generalized energy E , depending on the mixture parameters ρ, c 1 , α 1 , w k , s and A i k and we discuss its definition below. As soon as the generalized energy is defined, one can compute the thermodynamic forces ∂E (1a), (1d), (1e), (1g) and source terms in (1b), (1e), (1f), (1g): Parameters θ, θ 1 , and θ 2 characterize the rate of the mixture shear stress relaxation, and the rates of relative velocity relaxation and pressure relaxation respectively and can depend on state variables. Thus, the definition of generalized energy E gives us a governing PDEs for the solid-fluid mixture written in terms of phase parameters. The pressure p, the shear stress tensor σ i j and the temperature T are also computed via the generalized energy as: Note that the solution to system (1) satisfies the energy conservation law in accordance with the first law of thermodynamics. Also note that in numerical simulations, it is necessary to use the energy conservation law (4) instead of entropy balance law (1g). As stated in [23], system (1), (4) belongs to the class of SHTC (Symmetric Hyperbolic Thermodynamically Compatible) systems and describes dynamic processes in the solid-fluid mixture at finite-strains.
Our goal is to apply governing equations (1) to the description of processes in the deformed medium saturated by the compressible (viscous or inviscid) fluid. We shall do it assuming that the fluid volume fraction α 1 is identical to the porosity φ of the medium, α 1 = φ. Before discussing the definition of the generalized energy it is important to emphasize that we take the phase mass densities and mixture distortion as parameters characterizing the deformation of the element of the medium. This choice of parameters is made due to the fact that the parameters of the medium associated with volumetric deformation, such as density, energy and pressure, are additive values with respect to mass or volume fractions. As for the general strain, one can determine the strain tensor of each phase separately, but it is unclear how to formulate the equations for these individual strains and how to take into account their interaction. That is why we use the distortion of the entire mixture as the strain measure of the medium. Its evolution is determined by the velocity of the mixture via equation (1b).
To define the generalized energy density ρE , we, as in [23], assume that it is a sum of the kinetic energy ρE 0 of the center of mass of the mixture element, the kinematic energy of the relative motion ρE 1 , the energy of volumetric deformation ρE 2 , and the energy of shear deformations ρE 3 : The kinetic energy of the center of mass of the mixture element and kinematic energy of the relative motion are defined as Not that due to the definition of v j , w j , The energy of volumetric deformation is defined as The energy of shear deformation depends on the distortion of the entire solid-fluid mixture element and we define it as where c s,m is the shear sound speed in the mixture to be determined, and g is the normalized Finger (or metric) strain tensor: We assume now that the saturating fluid is viscous and its rheology is described by the unified SHTC (Symmetric Hyperbolic Thermodynamically Compatible) model of continuum [17,6]. This unified model considers a fluid as a visco-elastic medium with a small shear strain relaxation time. In such an approach, the solution to the Navier-Stokes equations of Newtonian fluids is approximated well by the solutions to the unified solid-fluid model [6] in the relaxation limit. Later, we will see that in the limiting cases of a pure fluid (α 1 = 1, α 2 = 0) or a pure solid (α 1 = 0, α 2 = 1), governing equations of the two-phase solid-fluid mixture (1) reduce to the unified SHTC governing equations for pure fluid and pure solid.
Because (9) is the shear energy of the mixture, the shear mixture sound speed c s,m should depend on the porosity (fluid volume fraction). We assume that c s,m is computed by the simple mixture rule via the shear sound speeds of fluid and solid phases c s,1 and c s,2 : Note that in [23], the two-phase model for deformed saturated porous medium was formulated under the assumption µ 1 = 0 that means that the saturating fluid is inviscid and has no resistance to shear.
Using the definition of the generalized energy (5), (6), (8), (9) and relations between mixture's and phase parameters', one can find thermodynamic forces E α , E ρ , E A k j , E w k , E c 1 and other thermodynamic parameters To close the model completely, one needs to define the parameters θ, θ 1 , θ 2 characterizing the rates of the mixture shear stress relaxation, relative velocity relaxation and pressure relaxation.
For further considerations it is convenient to formulate governing equations in terms of phase parameters of state, which can be obtained with the use of relationships between the mixture and phase parameters and formulae for thermodynamic forces (11)

Simplified isentropic single pressure model
In this section, we formulate a simplified model which is applicable to studying processes in porous medium with small temperature variations and instantaneous relaxation of interphasial pressures. If the saturated porous medium is initially in a thermal equilibrium, its deformations are small and there is no intensive thermal sources, then we can assume that the temperature variations in the process under consideration are small. This means that we can neglect thermal processes and can use the isentropic version of equations (12), neglecting equation (12g) and assuming that the phase energies and pressures do not depend on entropy.
Furthermore, we assume that the characteristic pore space in the deformed porous medium is small. This assumption allows us to reduce the set of state parameters and simplify system (12). Indeed, in this case we can assume that the pressure relaxation is instantaneous. The reason for this is that the phase pressures are equalizing due to the pressure wave propagating in the fluid in pores and skeleton. The pressure equalization time can be estimated as the time of several runs of the pressure wave in a pore of a characteristic size, and it is obvious that it is small in comparison with the characteristic time of interest to us in macroscopic samples. That is why we can assume that in model (12) the fluid and solid pressures are equal, p 1 = p 2 . Note that this assumption allows us to use the single pressure p = p 1 = p 2 in (12), and thus it is necessary to replace the equation for the volume fraction (12f) by the algebraic equation p 1 (ρ 1 ) = p 2 (ρ 2 ). We have to do this replacement because the relation p 1 (ρ 1 ) = p 2 (ρ 2 ) can be obtained as the relaxation limit of the model if θ 1 → 0. Taking into account the algebraic equation p 1 = p 2 , and using phase mass conservation equations (12c), (12d) one can derive the equation for the volume fraction of the single-pressure model: where K 1 = ρ 1 c 2 1 and K 2 = ρ 2 c 2 2 are the phase bulk moduli and c 1 and c 2 are the phase bulk sound velocities ( Thus, the simplified system takes the following form: The derivation of the governing equations for small amplitude wave propagation will be based on the above system, but first we note that in the limiting cases α 1 = 0, α 2 = 1 (pure solid) and α 1 = 1, α 2 = 0 (pure fluid) this system reduces to the unified solid-fluid model of finite-strain elastoplastic solid and viscous fluids respectively [17,6]. In fact, if we put α 1 = 0, α 2 = 1 in (14) and recall that v j = c 1 v j 1 + c 2 v j 2 , we then obtain the following system g i k g k j − 1 3 g l k g kl δ i j is the shear stress in the solid. System (15) does not include the limiting case of mass conservation equation for fluid phase (14d) because of the absence of the fluid (α 1 = 0). We see also that equations (15d), (15e) have no influence on the solid behavior, which is fully described by equations (15a)-(15c) because (15e) states the constancy of the phase volume fractions (α 1 = 0, α 2 = 1) and the relative velocity is not presented in all other equations. Thus the reduced model for the solid phase is governed by equations (15a)-(15c) which are exactly the equations of elastoplastic isentropic medium written in terms of velocities, pressure and shear stress, the full version of which, taking into account entropy, can be found in [6].
The similar system can be obtained in the limiting case of pure fluid. Thus, assuming α 1 = 1, α 2 = 0 and neglecting again the equation for the relative velocity, equation for the mass conservation of the solid phase, and for the volume fraction of the solid phase, we arrive to the following system (16) is similar to the isentropic version of the unified model of continuum presented in [17,6] and describes Newtonian viscous fluid flows as an asymptotic limit for small shear strain relaxation time τ if to choose that for small τ gives us This formula is exactly the definition of the Navier-Stokes shear stress and can be written as ∂x i is the strain rate tensor and η = ρτc 2 s,1 is the effective dynamic viscosity. The details of the asymptotic analysis and derivation of the Navier-Stokes stress for small τ can be found in Appendix A.
Thus, the presented hyperbolic thermodynamically compatible model describes a two-phase porous medium with a viscoelastic skeleton and a viscoelastic saturating fluid. This model allows us to derive differential equations for the propagation of small amplitude waves in a porous medium saturated with a viscoelastic fluid in the next section.

Governing equations for small amplitude wave propagation in an elastic porous medium saturated with a viscous fluid
In this section, a system of partial differential equations is derived that can be used to model small-amplitude waves in an initially unstressed medium. The derivation can be done by the standard linearization procedure and is similar to that presented in [23]. Consider a stationary unstressed medium and denote the parameters of state of this medium by the symbol "0". Assume that the initial value of the fluid and solid volume fractions α 0 1 , α 0 2 , (α 0 1 +α 0 2 = 1) are known.The immovability of the medium means that the mixture velocity and relative velocity in this initial state are equal to zero v i 0 = 0, w k 0 = 0 and therefore individual phase velocities are also equal to zero v i 10 = 0, v i 20 = 0. Furthermore, in the unstressed state, the pressure and shear stress of the mixture are also equal to zero, that means that p 0 Zero values of phase pressures correspond to the reference density values ρ 10 , ρ 20 and zero values of shear stress correspond to the distortion A 0 i j = δ i j . Now, assume that the solution has the form v k where the small perturbations of corresponding state variables are denoted by the symbol ∆. If to substitute (17) into system (14) and then neglect second and higher-order terms with ∆, we arrive to the following system of linear PDEs (c 0 s,m ) 2 = c 0 1 c 2 s,1 + c 0 2 c 2 s,2 , and ε i j = −(∆A i j + ∆A j i )/2 is the small deformation tensor. System (18) can be simplified by introducing new variables and reducing the number of equations. Let us take as a set of independent state variables the mixture velocity U i = (α 0 The coefficient θ characterizing the rate of stress relaxation is taken as θ = θ 0 τ, where θ 0 = 2(c 0 s,m ) 2 /ρ 0 , and τ is the relaxation time. The complete PDE system obtained from (18) and written in terms of the listed above variables reads as Here K is the bulk modulus of the mixture defined as µ is the shear modulus of the mixture which is defined according to (10) Note that in the complete system we do not need to include equation (18f) for the volume fraction (porosity) because its perturbation does not affect wavefields and, if necessary, ∆α 1 can be computed by solving this equation.
Thus, system of linear partial differential equations (19) has been formulated and it is applicable to modeling the propagation of small-amplitude waves in a stationary unstressed porous medium saturated with a fluid. It covers wave processes for the entire porosity range φ = α 0 1 ∈ [0, 1]. For the modeling of real media, it is necessary to determine the material constants in (19). In addition to the initially specified densities ρ 0 1 , ρ 0 1 and porosity φ = α 0 1 , we should define the bulk and shear modulus of the solid and fluid constituents K 0 1 , K 0 2 and µ 1 , µ 2 , which give us the mixture bulk K and shear modulus µ defined by (20), (21). It is also necessary to define coefficients θ 2 responsible for the interface friction and the shear stress relaxation time τ. In our study we take the the interfacial friction coefficient θ 2 constant, as in [23]. What concerns the shear stress relaxation time τ, we take it in the following form It is clear that the relaxation time should be a function of the porosity φ = α 0 1 , and this functional dependence is the subject of further theoretical studies. We take the empirical form of the relaxation time (22)  the nonlinearity on φ gives us a reasonable values of τ in the vicinity of pure phases. If the skeleton is pure elastic, then we take τ 2 = ∞ and the relaxation time τ is computed from the formula Fluid shear modulus µ 1 , relaxation time τ 1 and dynamic viscosity η 1 are connected by the relation Table 1 shows some values of the relaxation time (23) for several values of n and shear sound speed c s,1 of the fluid. All of the above-mentioned material constants should be chosen in such a way that the behavior of the medium under consideration corresponds to the available experimental data or the results known in the literature. In [23], a comparison of the two-phase thermodynamically compatible model of a porous medium saturated with an inviscid fluid and Biot's model is carried out. It is shown that the appropriate choice of material constants in the two-phase model allows one to obtain velocities of the fast and slow longitudinal waves and their frequency dependencies close to those obtained by Biot's model. To our knowledge, there are only a few research works in which the influence of the viscosity of a saturating liquid on wave fields in saturated porous medium is studied [26,8]. The main goal of these studies is the analysis of shear wave propagation and attenuation. In [26], it is discovered that the accounting of viscosity of the saturating liquid generates the so-called slow shear wave. Indeed, such a wave appears instead of a typical elastic shear wave and looks like a certain diffuse wave, the smeared structure of which arises due to the viscosity of the saturating liquid. In the above-mentioned papers, it is reported that a slow shear wave is rather difficult to observe, and it can be seen when elastic waves interact with interfaces between media with different porosities. In Section 7, shear waves in a porous medium saturated with a viscous fluid can be seen in some test cases, despite the strong attenuation caused by viscosity. As far as we know, there is no available experimental data concerning the shear wave propagation that is why we present qualitative study of such kind of waves. Since our goal is to see the main features of shear waves in a porous medium, we will choose material constants in such a way that clearly demonstrates these features. Then, as soon the experimental data will be available, one can vary material constants in order to achieve an agreement between theoretical and experiments studies.
The roots of the characteristic polynomial det(A − λI) are given by the formulas The speeds λ 1,2,3,4 correspond to the left-propagating and right-propagating so-called fast C fast and slow C slow characteristic speeds of compression waves (P-waves), while λ 5,6,7,8 are the shear wave characteristic speeds C shear . Fig. 1 shows the dependence of C fast , C slow , and C shear on the porosity φ = α 1 .

Dispersion relations
If the dissipation mechanisms (fluid viscosity and interfacial friction) are presented in (19), i.e. if the relaxation parameters τ < ∞ or θ 2 < ∞ are finite, then propagating waves exhibit dispersion, that is the speed of a wave (including small amplitude waves) depends on the wave frequency ω. For example, as shown in [15], the speeds of small amplitude waves V fast (ω), V slow (ω), and V shear (ω) in a medium with relaxation processes are smaller than the characteristic speeds (28) in the same medium but without dissipation mechanisms. In the high-frequency limit (ω → ∞), the sound speeds tend towards the characteristic speeds (28).
One can notice the following differences between two cases presented in Fig. 2 and Fig. 3. First, the shape of the phase velocity V slow (ω) changes from a longitudinal-like shape in Fig. 2 to a shear-like shape in Fig. 3. The low frequency limits (ω → 0) of V slow (ω) are different as well, 509.2 m/s in Fig. 3 versus 0 m/s in Fig. 2. Second, the low frequency limits of V fast (ω) are also very different. Third, the attenuation factor a λ (ω) of the fast mode increases over 3 orders of magnitude from Fig. 2 to Fig. 3. And finally, the attenuation factor a λ (ω) of the slow mode increases even more dramatically over 6 orders of magnitude from Fig. 2 to Fig. 3 Table 2.

General descriptions
In this section, we numerically verify the correctness of the obtained theoretical results and compare computed wave velocities for different frequencies with the velocities corresponds to dispersion curves. To generate seismic waves, it is necessary to specify the source term in the right-hand side of system (19). On can define a volumetric-type source exciting only compressional waves in a homogeneous media by inserting source term in the pressure equation (19c) or in equations for the diagonal components Σ i i (i = 1, 2, 3) of the deviatoric stress (19d). To generate all types of waves propagating from the source in a homogeneous medium, including S wave, it is necessary to define oriented-type source by adding source term only in one of the equations for Σ i i . In this case, we obtain vertical-or horizontal-type sources with maximum seismic energy propagation in the corresponding direction. We define the source term as where δ-Dirac's delta function, localizing the source in point (x 0 , y 0 ) and f (t )-Ricker's wavelet in time where ω = 2π f 0 is the angular frequency, f 0 is peak frequency in hertz and t 0 is the time wavelet delay, selected as t 0 = 2/ f 0 in our further consideration.
As a computational method, we choose finite difference schemes on staggered grids. The equations (19) form a symmetric hyperbolic system which is well suited to using staggered grids. The method is widely used due to computational efficiency and easy implementation. We use the classic second-order accurate staggered-grid stencil designed on the Cartesian coordinate system. The derivation of finite-difference approximation for system (19) can be found in [23]. To suppress artificial reflections from non-physical computational boundaries, we apply the perfectly matched layer technique (PML) in the original split-field formulation [5].

Homogeneous medium
We consider a set of two-dimensional test cases, in which the wave propagation is studied in a square domain Ω covered by a numerical grid with N grid points in each spatial direction and grid space interval ∆l so that Ω = [N ∆l ] 2 with centering at the origin of the Cartesian coordinate system (X , Y ). To simulate an unbounded medium, the domain Ω is surrounded by PML layers to avoid non-physical reflections from computational boundaries.
We consider three typical angular frequencies ω = 10 5 , 10 3 and 10 [rad/s], covering all cases of the characteristic behaviour of the dispersion curves presented in the previous section. For each selected angular frequency ω, we scale the spatial step ∆l and recording time T by some scaling parameter ε, leaving the number of nodes N constant. This is done because the seismic wavelength is proportional to the frequency and hence the size of the computational domain and the recording time must change with the frequency scaling. The parameter ε is selected for the convenience of comparison: if the wavefields properties do not depend on the frequency, then for all three cases the seismograms should be identical in scales proportional to this parameter.
First, let us consider the case ω = 10 5 for which the velocity dispersion curve reaches the upper asymptotic constant. For simulation we take grid parameters N = 10 4 , ∆l = ∆l 0 = 5 · 10 −4 m, medium parameters from Table 2, porosity φ = α 0 1 = 0.2, power index n = 8, vertically-type source location (x 0 , y 0 ) = (0, 0) and seismogram recording time T = T 0 = 10 −3 s. Grid spacing ∆l 0 provides approximately an amount of 10 points per slow compressional wavelength that is needed to avoid the numerical dispersion. Fig. 4 (left) shows a snapshot of the normalized total velocity vector for time t = 5 · 10 −4 s. The appearance of fast and slow compressional P waves and shear S wave exited by a vertically-type source radiating in the Y -direction can be observed. It is clearly can be seen the difference in radiation patterns, denoting the angular dependence of the strength of the waves radiated by a source for P and S waves: the maximum amplitude of P wave corresponds to the Y -direction, while for S wave this direction is ±45 • . Fig. 4 (right) plot the radiation patterns for the P (red line) and the S (green line) waves. The amplitude of the S wave is totally vanishing in directions that coincide with the X and Y grid axis. This remark explains why we choose the diagonal of the computational domain as a seismogram recording line. For all subsequent tests in this section, we assume that the receivers locate on the diagonal of Ω with a uniform spacing (dotted white line in Fig. 4).
For each selected angular frequency ω = 10 5 , 10 3 , 10 we simulate wavefields with different values of interfacial friction (relative velocity relaxation) coefficient and shear stress relaxation time to study their influence on the wavefield formation. In case of ω = 10 3 we set T = εT 0 s and ∆l = ε∆l 0 m with ε = 10 2 and for the case of ω = 10, we set ε = 10 4 . Fig. 5-7 show seismograms for the mixture vertical velocity v 2 for different ω and different values of relaxation parameters: relative velocity relaxation time θ 2 and shear stress relaxation time τ. First we examine the case when the medium does not have any dissipation (τ = θ 2 = ∞). Then we consider the case with two dissipation mechanisms having the same finite values of shear stress and relative velocity relaxation times (τ = θ 2 = 3.75 · 10 −4 ). The value of shear stress relaxation time is obtained by formula (24) for the viscosity η 1 = 10 −2 Pa·s and parameters n = 8, c s,1 = 100 m/s (Table 1). Additionally, in order to study the effect of the relative velocity relaxation only, we examine cases with τ = ∞ and two different parameters of θ 2 .
The comparison of seismograms in Fig. 5-7 allows us to formulate a conclusion about the dependence of wavefields on the two dissipation mechanisms. It is clearly seen that in the absence of dissipation, the velocities of the compressional and shear waves does not change for all frequencies. The situation is different when we switch on both dissipation mechanisms (case τ = θ 2 = 3.75 · 10 −4 ). The dependence of velocity on frequency that coincides with the theoretical one predicted by the dispersion curves in Fig. 1 can be observed. In addition, for frequency ω = 10 3 we see rather strong wave attenuation due to the predicted attenuation factor depicted in Fig. 2. Comparing seismograms for all three cases (τ = θ 2 = 3.75·10 −4 ), (τ = ∞, θ 2 = 3.75·10 −4 ) and (τ = ∞, θ 2 = 3.36·10 −7 ) we can see that the wavefield attenuation mostly depends on parameter τ and the parameter θ 2 is responsible for the appearance and attenuation of the Biot mode.

Two-layered medium
In this section, we consider a two-layered elastic/poroelastic medium with a horizontal interface located in the origin x = 0. The purpose of this experiment is to study the manifestation of a disperse shear S wave near the layer interface. Let us suppose that the upper part of the computational domain Ω is a pure elastic medium (φ = 0) with physical parameters Solid1 from Table 3 whereas the lower part of Ω is a poroelastic medium with the porosity φ = α 0 1 = 0.2 and material parameters Solid2, Viscous Fluid2 from Table 3. The size and numerical grid characteristics of Ω are the same as in the previous homogeneous case for ω = 10 3 . The volumetric-type source with frequency ω = 10 3 H z is located in the upper half-plane at point (x 0 , y 0 ) = (0, −30) m and it excites only one compressional wave in an elastic medium. The choice of volumetric-type source is explained by the desire to avoid the complex wavefield mixture and to focus on the analysis of reflected and transmitted waves generated by only one incident P-wave. For comparison we consider two cases of poroelastic layer with τ = θ 2 = ∞ and τ = θ 2 = 3.75 · 10 −4 . Fig. 8 and Fig. 9 present the results of numerical simulations. Fig. 8 shows wavefield snapshots of the normalized total velocity vector at time t = 8 · 10 −2 s and Fig. 9 shows seismograms for time T = 10 −3 s recorded in uniform spacing receivers located on the vertical line shifted 130 m right from the source position (dotted white line in Fig. 8). It clearly can be seen the wavefield difference for this two tests. Although the parameters of the upper elastic layer do not change, the wavefield in the upper layer changes due to reflection from the boundary of the poroelastic layer with different physical characteristics. If to take into account the shear stress and velocity relaxation times τ = θ 2 = 3.75 · 10 −4 , we can observe the occurrence of the disperse shear S wave with low velocity in the poroelastic layer. The amplitude of this wave attenuates rather quickly, so we can observe it near the interface only, where the wave appears and then attenuates. The shear wave can be seen more clearly if the amplitude of the wavefield is increased by two orders (see Fig. 8 1500 100 1040 10 −2 τ 1 = η/µ 1 8 Table 3.: Material parameters used in the layered model .

Model validation and verification
Note that the formulated model (19) can be used for numerical simulations of wave propagation for a quite wide range of problems. It can serve for wavefields simulation in a porous medium saturated by viscous and inviscid fluid for all range of porosity φ = α 0 1 ∈ [0, 1]. In this section, we consider test problems illustrating the applicability of the model for wavefield simulation in the medium containing regions of pure isotropic elastic solid φ = 0 and pure fluid φ = 1 with the single PDE system (19). We compare numerical results obtained by solving equations (19) with those obtained by the classical approach based on the stress-velocity elasticity model with interface tracking between elastic medium and fluid, realized by the boundary conditions treatment in the numerical method. Numerical examples presented in this section include wavefields simulation in a finely layered medium composed of elastic isotropic solid/solid or solid/fluid alternate layers.
We study a periodic, two-layered system inside a square computational domain [500m] × [500m] centered at the origin with 5000 2 grid mesh and uniform grid spacing 0.1m. Layers of the thickness 0.5m are settled horizontally. For linear elasticity model each layer is defined by a prescribed boundaries and the medium inside is defined by the density ρ i (i = 1, 2), the compressional velocity c p,i (i = 1, 2) and the shear velocity c s,i (i = 1, 2). For two-phase model (19) the interface between layers is defined via the porosity function φ. In the absence of the velosity and shear stress relaxation (τ = θ 2 = ∞), φ = 0 corresponds to the pure elastic medium and φ = 1 corresponds to the pure fluid.
We consider two cases of finely layered medium with parameters Solid1/Fluid1 and Solid1/Solid3 from Table 3. Fig. 10 shows snapshot of the total velocity vector computed for the Solid1/Fluid1 layered medium at time t = 6.2 · 10 −2 s. The volumetric-type source with dominant frequency f 0 = ω/2π = 10 3 /2π H z is located in a fluid layer at the center of the domain. For selected frequency and medium parameters, the thickness of the layers corresponds to 1/75 of the dominant P-wavelength in the solid. To illustrate the agreement between the solutions obtained by the linear elasticity model and the two-phase model (19), Fig. 10 is split to show both numerical results. On the left side, the snapshot obtained by the two-phase model (19) is depicted, whereas on the right side the snapshot obtained by the classical stress-velocity formulation of the linear elasticity equations is shown. It can be seen that they are practically identical, which is confirmed by semblance value 99,9% computed by the formula where a i and b i are the wavefield components changing with time in both simulations. Good agreement of the results is achieved by a practically identical finite-difference schemes approximating the selected model for both calculations.
The similar excellent agreement between solutions obtained by the linear elasticity model with interface tracking and the two-phase model can be obtained if to consider the case of elastic finely layered medium with layers parameters Solid1/Solid3. In this case for wavefields simulation we just replace Fluid1 layers in the previous consideration with the Solid3 layers. Fig. 11 (a) shows a two-piece snapshot of the total velocity vector computed by the model (19) and linear elasticity equations for the layered medium Solid1/Solid3 for the source frequency ω = 10 3 at time t = 6.2 · 10 −2 s.
A significant difference of wavefield in Fig. 10 and Fig. 11 can be explained by the homogenization theory for a finely layered medium composed of periodic elastic layers. The stack of the horizontally parallel thin isotropic layers acts as a VTI (transverse isotropy with a vertical axis of symmetry) anisotropic medium in the long-wavelength limit. We use Backus averaging procedure [1] to estimate the so-called quasi-P wave and quasi-S wave phase velocities in order to compare with our numerical results. To this end, the position of the peak amplitude of quasi-P (white dashed curve) and quasi-S (yellow dashed curve) wavefronts is imposed on the wavefield snapshot in Fig. 11(right). We remind that in the numerical modeling, the source wavelet is shifted in time by 2/ f 0 = 1.2 · 10 −2 s, that is why the phase velocity curves are calculated at t = 5 · 10 −2 s after the source excitation. We can observe that these curves are in good consistency with the finite-difference results, thus demonstrating the correctness and accuracy of our model.
At the end of this section, we present one more calculation for a finely layered medium containing poroelastic layers with a viscous fluid in order to investigate the influence of the fluid viscosity on the wavefield behavior. For comparison purposes, we slightly modify the previous experiment with elastic layers Solid1/Solid3, and substitute the pure elastic layer Solid3 with a poroelastic layer Solid3/Viscous Fluid2 from Table 3 with porosity φ = 0.2. All other parameters are the same. Fig. 12 shows snapshot of the normalized total velocity vector computed for this medium at the times t = 1.5 · 10 −2 s and t = 6.2 · 10 −2 s. Comparing Fig. 11 and Fig. 12 demonstrates a significant change in wavefield behaviour produced by adding porosity and viscosity into the model. We observe strong seismic attenuation of all waves in vertical direction caused by the combined effect of multiple inter-layer reflection/refraction scattering and viscous mechanism. Also, note that wavefield front doesn't gather in the horizontal direction.

Strongly heterogeneous medium
This experiment is carried out for the porous media model downloaded from Comsol Application Gallery, Pore-Scale Flow, produced by scanning the electron microscope images during the pore-scale flow experiments and schematically presented in Figure 13. We use this geometry to define porous medium with the material parameters from Table  2 and porosity φ = 0.8, containing less permeable porous channels with porosity φ = 0.2 ( white color in Fig. 13). For numerical simulations let us consider two cases of relaxation mechanisms with τ = θ 2 = ∞ and τ = θ 2 = 3.75 · 10 −4 . Fig. 14 show wavefield snapshots of the total velocity vector computed for ω = 10 4 at times t = 1.2 · 10 −3 s and t = 2.4 · 10 −3 s. At time t = 1.2 · 10 −3 s, we see the propagation of the P wave and the initiation of the Bio wave for both cases. But in the case of τ = θ 2 = 3.75 · 10 −4 , the wavefield's amplitude is strongly attenuates from the beginning and makes it impossible to form the Biot mode ( see time t = 2.4·10 −3 s ) if the relative velocity and shear stress relaxations are taken into account.

Conclusions
We have presented a new hyperbolic two-phase model of a porous medium saturated by a viscous fluid. The governing equations of the model are derived by generalizing the unified model of continuum [17,6], which in turn is derived using the theory of Symmetric Hyperbolic Thermodynamically Compatible (SHTC) systems [11,21,16]. The unified model of continuum is based on the first principles, mathematically well-posed, satisfies thermodynamics laws and describes heat-conducting solid and liquid states of a medium by a single hyperbolic PDE system. Consideration of the saturated porous medium as a two-phase solid-fluid mixture and coupling of the unified model with the twophase SHTC model of compressible fluid [19] allows us to design the SHTC nonlinear model for compressible fluid flow in a deformed porous medium. The presented model, similarly to the unified model, can take into account the viscosity of a saturating fluid in a hyperbolic formulation. Thus, in addition to interfacial friction, there is a dissipation in the model caused by the viscosity of the saturating fluid which should be taken into account in the interpretation of seismic data. Since there are no available experimental or well-established theoretical data on the properties of shear waves in the literature, we choose some empirical functions of state parameters for closing relations such as the equation of state (specific internal energy), interfacial friction coefficient and shear stress relaxation function.
Using the presented nonlinear finite-strain SHTC model, the governing equations for the propagation of smallamplitude waves in a porous medium saturated with a viscous fluid are derived. As in the well-known theories of porous media, three types of waves can be found: fast and slow compression waves and shear wave. It turns out that the shear wave attenuates rapidly due to the viscosity of the saturating fluid, and this wave is difficult to see in typical test cases. However, some test cases are presented in which shear waves can be observed in the vicinity of interfaces between regions with different porosity.
The presented model can be validated by comparison with experimental data on the propagation of waves in porous media and, in particular, on the propagation of shear waves, or by modeling flows in an elastic skeleton on a pore scale, which is one of the goals of our forthcoming studies.
From (43) one can see that the zero term G 2 0 should satisfy the uniform volumetric deformation, i.e. G 2 0 − 1 3 tr(G 2 0 )I = G 2 0 −ΛI = 0, where Λ = (det G 0 ) 2/3 (the latter follows from the chain of equalities det G 2 0 = det(G 0 ) 2 = Λ 3 ). Note that the density ρ 1 is a product of the reference density ρ 10 and square root of det G: ρ 1 = ρ 10 det G, but the two first terms of its expansion (zero and first ones) has no influence on the final result and we use the notation ρ 1 everywhere below. Now, with the use of equation (43) we can derive equation for G 2 : Substituting (42) into (43), we obtain for the first order terms in ε Then, taking into account that for zero terms G 0 = Λ 1/2 I and Substituting the above formula and G 0 = Λ 1/2 I in (44), we obtain Finally, for the first term of the asymptotic expansion, we obtain where˙ i k = ∇v 1 + ∇v T 1 is the strain rate tensor and η = τc 2 s,1 is the dynamic viscosity.     Snapshot of the normalized total velocity vector for ω = 10 3 at time t = 8 · 10 −2 s computed with parameters τ = θ 2 = ∞ (left column) and τ = θ 2 = 3.75 · 10 −4 (right column) for two-layered medium. Notations: black line -layer interface, white dotted line -receivers location, 1 = direct P-wave generated by the source, 2 = reflected P-wave, 3 = reflected S-wave converted from the P-wave, 4 = P-to-S head wave, 5 = transmitted fast P-wave, 6 = converted and transmitted S-wave, 7 = converted and transmitted slow P-wave.