Wave Equations of Porous Media Saturated With Two Immiscible Fluids Based on the Volume Averaging Method

Simulating and predicting wave propagation in porous media saturated with two fluids is an important issue in geophysical exploration studies. In this work, wave propagation in porous media with specified structures saturated with two immiscible fluids was studied, and the main objective was to establish a wave equation system with a relatively simple structure. The wave equations derived by Tuncay and Corapcioglu were analyzed first. It was found that the coefficient matrix of the equations tends to be singular due to the inclusion of a small parameter that characterizes the effect of capillary stiffening. Therefore, the previously established model consisting of three governing equations may be unstable under natural conditions. An improved model based on Tuncay and Corapcioglu’s work was proposed to ensure the nonsingularity of the coefficient matrix. By introducing an assumption in which one fluid was completely wrapped by the other, the governing equation of the wrapped fluid was degenerated. In this way, the coefficient matrix of wave equations became nonsingular. The dispersion and attenuation prediction resulting from the new model was compared with that of the original model. Numerical examples show that although the improved model consists of only two governing equations, it can obtain a result similar to that of the original model for the case of a porous medium containing gas and water, which simplifies the complexity of the calculations. However, in a porous medium with oil and water, the predictions of dispersion and attenuation produced by the original model obviously deviate from the normal trend. In contrast, the results of the improved model exhibit the correct trend with a smooth curve. This phenomenon shows the stability of the improved model and it could be used to describe wave propagation dispersions and attenuations of porous media containing two immiscible fluids in practical cases.


INTRODUCTION
Wave propagation in partially saturated porous media is of interest in the research areas of geophysics, petroleum engineering and underground science. For instance, partial saturation of two or more fluid components usually occurs in complex unconventional reservoirs (Santos et al., 2019).
There are plenty of works related to this issue. Some scholars have studied modulus calculations for composite materials. White (1975) established the first bulk modulus expression that considered dissipation mechanisms on a mesoscopic scale, which is also known as the patch saturation model. A model of concentric spheres with periodic distribution was formulated to simulate the heterogeneity of fluids in porous media. Then, based on the quasistatic Biot theory (Biot, 1956;Biot, 1962), Johnson (2001) proposed a dynamic frequency-dependent bulk modulus to calculate the dispersion and attenuation relations of patches with arbitrary geometries. The bifurcation function of the bulk modulus under different limit conditions was considered. Subsequently, Tserkovnyak and Johnson (2003) improved this model by considering the effect of membrane surface tensions between two fluids in pore space. Thereafter, a patch saturation model with a random distribution was proposed by Müller and Gurevich (2015).
In the work of wave equations, there are two main modeling methods: the mixture theory and volume averaging theory. It is assumed in the mixture theory that the pore media is homogeneously distributed on a macroscopic level and it is therefore not necessary to describe the pore structure (Tuncay and Corapcioglu, 1997). There is a vast body of work operating under the mixture theory, and only a few of them are listed here. In early studies, Brutsaert and Luthin (1964) extended Biot's theory (Biot, 1956;Biot, 1962) to a case with two fluids, and three kinds of compressional waves were predicted. Based on the principle of analytical mechanics, Berryman (1986), Berryman (1988) first established the kinetic and potential energy functions of a porous medium containing two fluids. Then, wave equations for partially saturated media based on separating fluid were established to analyze the dissipation caused by the interactions between fluids and solids. The model inadequately considered the influence of the interaction between the two fluids on the seismic wave energy . According to the fluid separation theory (Berryman, 1986;Berryman, 1988), as mentioned above, Santos et al. (1990a), Santos et al. (1990b), Santos et al. (2004) derived wave equations for porous media with two fluids. The capillary pressure function was introduced, and the idealized test method was used to determine the elastic coefficients. However, some parameters of the model may have had no practical physical meaning due to the use of the idealization test . Lo et al. (2005), Lo et al. (2007), Lo et al. (2015) further studied the model proposed by Berryman (1986), Berryman (1988) and added an inertial coupling term in the equations. Ba et al. (2016) studied the effect of rock matrix stiffening caused by clay squirt flow on wave dispersion. Moreover, Ba et al. (2017) proposed a double double-porosity model to describe wave propagation in anelasticity rock. Such a model takes into account both the heterogeneity of rock fabric and that of fluid distribution.
Another important way to derive wave equations is the volume averaging method, which benefits from the proof and development of the volume averaging theorem (Anderson and Jackson, 1967;Slattery, 1967). It does not require the assumption of a uniform spatial distribution of pores. The variables have clear physical meanings, and there is no need to consider fluid distribution shapes on a microscopic level with this technique. In this way, motion equations and the constitutive relations can be obtained by taking the average of the corresponding microscopic equations. A theorem for the volume averaging of a gradient is employed during the derivation process.
De la Cruz and Spanos (1983), De la Cruz and Spanos (1985) used the volume averaging method to derive wave equations for porous media saturated with one fluid. It was demonstrated by numerical examples that the results obtained from their final equations are almost the same as those of Biot's theory within a certain frequency range for the case of a given fluid with a relatively low viscosity. The volume averaging method was also used in the modeling of the so-called double-porosity dual-permeability model developed by Pride and Berryman (2003a), Pride and Berryman (2003b). They used this method to define the physical quantities in a heterogeneous porous media. A representative work of wave equations considering two fluids is a model developed by Tuncay and Corapcioglu (1995), Tuncay and Corapcioglu (1996), Tuncay and Corapcioglu (1997). They applied this method to explore wave propagation in poroelastic media saturated by two immiscible fluids. The capillary effect was considered, and they proposed a new way to express the relation of volume variation to determine the coefficients of the constitutive relations. Alternatively, one can use the two-space method proposed by Burridge and Keller (1981). They started from the microscopic, linearized elasticity equations for the solid phase and the Navier-Stokes equations for the fluid component. The asymptotic expansion technique was used. Although the derivation process was complicated, the final equations they obtained were equivalent to Biot's equations when the viscosity of the fluid was low. In principle, the final expressions derived by the volume averaging method are equivalent to those of the two-space method (Tuncay and Corapcioglu, 1997).
Comparing the two mainstream methods, the former assumes a porous medium with a homogeneous structure, regardless of the geometric details. The specific geometric structure of the pore space is considered by the latter. Based on this difference, we think the volume averaging method is a more general modeling method since it is convenient to obtain a physical model in line with the real world.
This study aims to establish numerical stable wave equations for porous media containing two fluids using the volume averaging method, and the manuscript is organized as follows. In Section 2, the model describing wave propagation in porous media saturated with two immiscible fluids developed by Tuncay and Corapcioglu is analyzed. It is found that the coefficient matrix of the equations tends to be singular due to the small parameter that characterizes the capillary stiffening effect. Then, wave equations for porous media saturated with an effective fluid are proposed as an improved model based on their study. Two numerical examples related to dispersion and attenuation are given in Section 4. From the results, the improved model can obtain similar results as the original model in the case of porous media containing gas and water. However, for the case of oil and water, the curve calculated by the original model is unstable; one can obtain smooth solutions using the improved model.

The Analysis of Tuncay and Corapcioglu's Model
Considering a porous medium containing two immiscible fluids, Tuncay and Corapcioglu (1997) started from a microscopic momentum balance equation between the solid and the two fluids, respectively. Only the low-frequency case was considered in their study. The corresponding macroscopic equations were obtained by the volume averaging method. Meanwhile, they obtained the macroscopic constitutive relation in a novel way (Tuncay and Corapcioglu, 1997). The final wave equations are obtained by combining the macroscopic momentum balance equation and the constitutive relation, which can be written as follows (Tuncay and Corapcioglu, 1997): Where u s is the solid displacement and u f 1 and u f 2 are the two fluid displacements. Note ϕ is the porosity and S 1 and S 2 are the saturations of the two fluids in which S 1 + S 2 1 and ϕ ϕ 1 + ϕ 2 , where ϕ 1 S 1 ϕ and ϕ 2 S 2 ϕ. ρ s , ρ f 1 and ρ f 2 are the mass densities of the solid and the two fluids, respectively. The averaged densities are denotes as 〈ρ s 〉 (1 − ϕ)ρ s , 〈ρ f 1 〉 ϕ 1 ρ f 1 and 〈ρ f 2 〉 ϕ 2 ρ f 2 . C 1 and C 2 are two dissipation coefficients expressed as C 1 ηϕ 2 1 /κ 1 and C 2 ηϕ 2 2 /κ 2 , where η denotes the fluid viscosity; and κ 1 and κ 2 are the relative permeabilities of the two fluids in pore space, respectively. Note that t represents the time, and the dot above the physical quantity indicates the derivative of time. In addition, the elastic constants are expressed as follows (Tuncay and Corapcioglu, 1997): A 1 , A 2 , and A 3 are expressed as follows: Where K s , K f 1 and K f 2 are the bulk moduli of the solid matrix material and two fluids, respectively; K b and G are the drained bulk modulus and shear modulus of the solid skeleton, respectively. Note a p 11 a 11 + G/3. β and n are two fitting parameters needed to calculate the capillary force p cap (Van Genuchten, 1980;Tuncay and Corapcioglu, 1997;Lo et al., 2005).
Equation 1 contains three governing equations, which is the major difference between the wave equations of porous media saturated with two fluids and that of one fluid. It is obvious that the structure of the second equation is the same as that of the third equation, but the coefficients in front of each term differ.
It is worth noting that an important contribution of this model was that the effect of the capillary force was considered through the elastic constant A 2 , which is generally a positive real number. As seen in Eq. 2, A 2 interacts with elastic parameters, such as K f1 , and characterizes the strengthening effect of the capillary force on the elastic parameters. However, as will be shown in Section 3, the magnitude of A 2 , as calculated by Eq. 3, is much smaller than that of K f1 and K f2 , especially in the case of a porous medium saturated with two liquids. Different works (Tserkovnyak and Johnson, 2003;Pride et al., 2004;Kobayashi and Mavko, 2016) have studied this issue. For example, Kobayashi and Mavko (2016) showed that the effect of capillary forces on the P-wave velocity was not significant by testing a set of experimental samples. Thus, it seems that the strengthening effect of the capillary force on elastic parameters can be ignored in most cases since A 2 has a small value compared with the related parameters and can be neglected.
Once the relations of A 2 ≪ K f 1 , A 2 ≪ K f 2 and A 2 ≪ K s are identified, one can represent the relationship between the elastic coefficients of the second and third governing equations as follows: Equation 4 shows that the three elastic coefficients of the second governing equation are proportional to that of the third equation.
In other words, the matrix composed of these nine coefficients tends to be singular. This may cause the instability of the solutions of Eq. 1 since the determinant of the matrix will appear in the Christoffel equation of the compressional waves (Tuncay and Corapcioglu, 1997). In this work, numerical examples are given to illustrate this phenomenon, and the corresponding mathematical proof of instability needs to be studied further. Generally, it is necessary to consider the influence of capillary forces in the study of wave propagation in porous media containing two fluids, which depends on the characteristics of the two fluids and the conditions. Therefore, the influence mechanisms of capillary forces on wave propagation need to be further studied, for example, to establish a new calculation method.

The Effective Fluid Model
The above analysis shows that this potential issue can be eliminated by combining the latter two governing equations. Hence, one solution is to reduce the number of equations from three to two by equivalent reduction of Tuncay and Corapcioglu's model. First, we review the derivation of the Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 618909 constitutive relation for a porous medium with two immiscible fluids. Three relations were provided as follows (Tuncay and Corapcioglu, 1997): 1) The definition of the capillary force is as follows: Where δS 1 represents the disturbance to S 1 .
2) The bulk deformation of the solid phase, which can be expressed as follows: Where p s , p f 1 and p f 2 are the pressures exerted on the solid and two fluids, respectively. The bar above each quantity indicates that it is obtained by averaging the volume of the microscopic quantity.
3) The volume variation relations of the representative element volume (REV) are calculated as follows: Where j s, f 1 and f 2 represents the relations for the solid and two fluids, respectively; α j is the volume fraction; P j is the effective stress of each phase; τ j is the stress tensor; K j is the bulk modulus; and u j is the macroscopic displacement of each phase. The expressions of the solid and fluid pressure can be solved according to Eqs. 5-7. The constitutive equation is obtained, and the expressions of the elastic coefficients are also determined at the same time (Tuncay and Corapcioglu, 1996;Tuncay and Corapcioglu, 1997).
Then, an attempt to transform the two fluids into an equivalent fluid in a reasonable way is made. Based on the results of Pride et al. (2004), a special case is considered in which one fluid phase (noted as fluid 2) is completely wrapped by the other phase (noted as fluid 1). Specifically, fluid 2 has no contact the solid skeleton. This means that there is no interaction between fluid 2 and the solid phase, while fluid 1 is in contact with both fluid 2 and the solid skeleton, which resembles the model constructed by White (1975). The discontinuous interface of fluid pressure still exists, and therefore, capillary forces still exist. Compared with the conventional effective fluid model, the main advantage of this effective method is that it retains the nonuniform fluid pressure distribution in space; as a result, the wave equation contains two independent fluid pressures.
Furthermore, the wrapped fluid (fluid 2) is assumed to have no contact with the boundary of the REV selected by the volume averaging method. Thus, the deformation of fluid 2 does not affect the whole unit, and one can obtain ∇ · u f 2 0. In this way, Eq. 7 with j f2, which represents the pressure of fluid 2, can be solved and expressed in terms of the related quantities. The resultant equation system can be listed as follows: The expressions of P s and P f 1 in terms of ∇ · u s and ∇ · u f 1 can be obtained by solving Eq. 8. Then, similar to Tuncay and Corapcioglu's definition (1997), one can obtain the following: The new elastic coefficients are denoted asã ij (i, j 1, 2) with the following expressions: where A 1 , A 2 and A 3 are the same as Eq. 3. Then, the constitutive relation of the effective fluid model is obtained as follows: where 〈σ〉 and 〈s〉 are the stress tensors of the solid and the effective fluid, respectively, and I is the unit tensor. Note that in Eq. 11, u f rather than u f 1 is used to denote the displacement of the effective fluid. Here, the effective fluid means that only fluid 1 is needed when the interaction between solid and fluid is considered, while the parameters of fluid 2 are needed for the calculation of the coefficients in the constitutive relation.
It is easy to verify that the elastic constantsã ij (i, j 1, 2) in the constitutive relation of the effective fluid model are equal to those of Eq. 2. In the following text, the elastic constants are uniformly noted as a ij (i, j 1, 2) for convenience.
The bars over u s and u f can be left out for simplicity. After combining Eqs. 11, 12, the final wave equations of porous media saturated with effective fluid are obtained as follows where C ηϕ 2 /κ. Compared with Eq. 1, the proposed model Eq. 13 is simpler in form since it includes only two governing equations. It is also easier to solve such equations. Although this model looks like that of the wave equations for a porous medium saturated with a single fluid, the expressions of the elastic coefficients are determined by two fluids. In addition, it is also similar to but different from Biot's equations (Biot, 1956;Biot, 1962) since there are no inertial terms of acceleration between the solid and fluid phases. Therefore, a relatively simple and stable wave equation structure is established, and a more complicated mechanism of wave propagation dissipation needs to be further proposed.

NUMERICAL EXAMPLES
Two sets of rock parameters are selected to verify the effectiveness and reliability of the improved model. The rock samples saturated with gas-water and oil-water are case 1 and case 2, respectively. Combined with the principle of plane wave analysis (Carcione, 2007), Eqs. 1, 13 are applied to calculate the attenuation and dispersion curves in these two cases. The detailed Christoffel equations for the compression and shear waves are given in the Appendix.

Porous Media Containing Gas and Water
A sample of Massillon sandstone (Tuncay and Corapcioglu, 1996) saturated with gas and water is selected for calculation. Here, fluid 1 is water and fluid 2 is air. The water saturation is denoted as S w . The experimental parameters are shown in Table 1. Four water saturations (S w 0.1, 0.5, 0.8, and 0.9) are chosen for comparison. Under the parameters of Table 1, the value of A 2 according to Eq. 3 is 326.56, which is two orders of magnitude smaller than K f2 and much smaller than K s and K f1 . Therefore, the strengthening effect of the capillary forces on the elastic parameters is not significant and can be ignored reasonably.
The dispersion curves of the two compression waves (fast and slow P-wave) and the shear wave (S-wave) are shown in Figure 1. For the dispersion curve of the fast P-wave, both sets of wave equations produce similar results at different S w values, as shown in Figures 1A,B. They both have close limits for high (>10 6 Hz) and low frequency (<10 4 Hz) range. Although each curve varies smoothly and presents an 'S' shape as a whole, it can be seen that the velocity curves resulting from the new model hold that in the intermediate frequency range; the greater the water saturation is, the lower the velocity. In contrast, there is no such consistency in the original model, and different characteristics occur in the intermediate-high frequency range. From previous work , it is known that the magnitude of velocity decreases with larger S w values in the low-frequency band, and the variations in velocity increase as S w increases. In addition, the velocity predictions of the fast P-wave by the Biot-Gassmann-Wood (BGW) model are 1,145.04 m/s, 1,150.34 m/s, 1,168.18 m/s, and 1,193.72 m/s, corresponding to the four saturations. They are basically consistent with the corresponding low-frequency limit values with different S w values predicted by the two models. The predictions by the Biot-Gassmann-Hill (BGH) model are 1933.24 m/s, 1769.30 m/ s, 1,461.24 m/s, and 1,239.29 m/s, corresponding to the four saturations, which are larger than the results predicted by the two models. This shows that the two models can obtain the velocities that are covered by two bounds; however, they still need to be improved further by introducing other physical mechanisms, such as local fluid flow.
For the dispersion curve of the slow P-wave, as shown in Figures 1C,D, one can see that the value of each curve is very small in the low-frequency band (<10 4 Hz) and then increases as the frequency increases, finally converging to a constant value in the high-frequency band. From Figure 1C, the higher the water saturation is, the larger the velocity will be, but they converge to the same value with increasing frequency. The results do not converge to the same value, as shown in Figure 1D. Different from the fast P-wave, the dispersion curves of the slow P-wave obtained by the two sets of equations are quite different compared with Figures 1A,B. This is because the number of governing equations is different, which leads to different mechanisms when slow P-wave propagation occurs in a porous medium. As we know that slow P-waves are not an elastic wave but rather a type of dissipative wave, its dispersion pattern differs from that of a fast P-wave. It should be noted that there is another type of slow P-wave for the solution of Eq. 1, while there is only one type of slow P-wave for Eq. 8, so the second type of slow P-wave is not shown here.  (Tuncay and Corapcioglu, 1996).

Parameter and unit
Symbol Value For the dispersion curve of the S-wave, as shown in Figures  1E,F, the results obtained by the two models are basically the same. The trend of the results is similar to that of the fast P-wave but with a lower magnitude.
According to the results of the fast P-wave and the S-wave under the sample condition, the results produced by the improved and original models are similar, especially within the frequency range of the seismic waves. This indicates that the original model developed by Tuncay and Corapcioglu can be replaced by the improved model proposed in this paper, and the improved model is simple in form and easier to solve. The velocity of the fast P-wave and S-wave vs. the water saturation were also calculated in their study (Tuncay and Corapcioglu, 1996), and the results matched well with the experimental data. We note that the predicted results of the improved model also match the experimental data well under the sample since one can obtain similar results as that of the original model.
The attenuation curves are shown in Figure 2. The attenuation curves of the fast P-wave and the S-wave resemble a downward parabola, and the maximum value is obtained at the frequency at which the velocity varies fastest in the corresponding dispersion curve. The attenuation curve of the fast P-wave is basically the same as that of the S-wave. In fact, the variation law of each attenuation curve can be roughly obtained from the dispersion curve. The two models have different trends in the low-and highfrequency bands. The frequency of the maximum attenuation is low for high S w values according to the original model, while for the new model, the opposite is true. As mentioned above, the velocity variations produced by the original model are irregular. Therefore, it can be considered that the results of the new model are more realistic.
Slow P-waves are usually a type of diffusion wave with a strong attenuation, whose attenuation curve is different from other wave types. It can be seen from Figures 2C,D that for each S w, the attenuation curve presents a horizontal line at low frequencies (<10 4 Hz). Then, this value gradually decreases linearly, and there is a smooth transition between the two lines. It should be noted from Figure 2C that the variation of the attenuation value for S w 0.9 is FIGURE 1 | The dispersion curve of seismic wave propagation in porous media containing gas and water: (A,C,E) are the dispersion curves of the fast P-wave, the slow P-wave and the S-wave corresponding to four saturations calculated by Eq. 1, respectively; (B,D,F) are those of dispersion curves by the improved model proposed in this work.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 618909 6 not a horizontal line like that of others in the intermediate frequency band (10 5 -10 7 Hz). Correspondingly, it can be seen from Figure 1C that the dispersion curve at S w 0.9 does not resemble the other curves since a convex shape appears in the curve, which is believed to be caused by the instability of the model. The instability of the model makes the solution unstable with variable parameters such as frequency.
In this example, the attenuation curves obtained by the two models are smooth overall, which is consistent with the results of some traditional models (Biot, 1956;Biot, 1962;Santos et al., 1990a;Santos et al., 1990b;Pride and Berryman, 2003a;Pride and Berryman, 2003b).

Porous Media Containing Oil and Water
The differences in the physical properties of oil and water are less notable than that of gas and water. To further verify the applicability of the established model, a rock sample containing oil and water is taken from a Columbian sandy loam unit (Johnson, 2001;Lo et al., 2005). The rock parameters are shown in Table 2. In this example, the magnitude of the intermediate parameter A 2 is 10.488, which is much smaller than that of the former case. Therefore, it is more reasonable to ignore the strengthening effect of capillary forces on elastic parameters for porous media containing two liquids according to Eq. 3.
The results of the dispersion curves are shown in Figure 3. It shows the dispersion curves of the fast P-wave calculated by the two models. It can be seen that the curves corresponding to different saturations no longer have an 'S' shape when using the model developed by Tuncay and Corapcioglu. The value remains constant within the low-and middle-frequency ranges (<10 6 Hz) at various water saturations, while the curves are convex at approximately 10 7 Hz and then decrease gradually. This trend is somewhat surprising and differs from that of case 1. As shown in Figure 3B, the curves obtained by the improved model remain 'S' shaped, similar to that in the case of gas and water. In this case, the model is more unstable due to the smaller A 2 value. Therefore, the mathematical solution of the unstable model can no longer. Represent the physical phenomenon correctly. It is difficult to obtain the velocity variation range of the results using the original model, and it is difficult to compare these results with those of the BGW and BGH models. In addition, the magnitude of the fast P-wave velocity in the former case is obviously lower than that resulting from the oil-water case as the physical properties of gas and oil are quite different, and the compressibility of oil is closer to that of water.
For the slow P-wave, as shown in Figures 3C,D, the trend of these results is similar to that of the previous example. However, it seems that the curve related to S w 0.5 in Figure 3C is not smooth enough. Note that the second type of slow P-wave in Eq. 1 is not listed here. Figures 3E,F show that the velocities of the S-waves corresponding to different saturations obtained by the two models are roughly the same. Differing from the fast P-wave,  (Johnson, 2001;Lo et al., 2005).

Parameter and unit
Symbol Value  Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 618909 8 the dispersion curves of the S-wave obtained by the two models are smooth. The corresponding curves are basically the same, but it can be seen from Figure 3F that the value tends to the same eventually while that is not the case in Figure 3E. This shows that the predictions output by the new model are close to those of the original model if its solutions are relatively stable.
The attenuation curve results are shown in Figure 4. It shows that the shape of the attenuation curve of the fast P-wave is very strange, which indicates that the model is unstable when predicting the attenuation of a porous medium containing oil and water. A smooth curve can still be obtained by the improved model, which suggests a stable solution. The attenuation value decreases with larger S w values. Figure 4C shows that the value of the attenuation curve is unstable in the high-frequency band. The attenuation curve of the S-wave is relatively stable, which can also be seen from the dispersion curve. These results indicate that the original model may have deficiencies when predicting the attenuation of oil and water.

CONCLUSION
Theoretical analysis showed that the model developed by Tuncay and Corapcioglu based on the volume averaging method may be unstable since the structure of the wave equations may be problematic. In this work, an effective fluid model of wave propagation is proposed based on their study. The dispersion and attenuation curves calculated by the former are not smooth enough in the case of a porous medium saturated with oil and water. In contrast, the effective fluid model can obtain smooth dispersion and attenuation curves, which are consistent with the results of classical models. The stability of the two models from a mathematical perspective as well as the corresponding numerical simulations need to be further studied. In addition, the improved model alone may not be enough to describe the mechanisms of wave attenuation and dissipation. It is still worth noting that wave equations with a simple and stable structure is found by introducing an effective fluid.
Future studies will focus on three points. Firstly, wave equations that consider more interaction mechanisms between the two fluids need to be developed. Secondly, it is necessary to compare the theoretical calculation results with more experimental data to verify the validity of the new model. Finally, the effect of the geometric structure of porous media on wave propagation should be studied. For example, a porous medium with porous space consisting of interconnected microtubes can be selected as the REV, and the permeability and elastic constants of such a network can be calculated (Xiong et al., 2020). Alternatively, one can consider cracked porous media with penny-shaped inclusions (Zhang et al., 2019) as the REV. Then, it is possible to obtain wave equations that connect the microscopic details of a porous medium with macroscopic wave propagation by the volume averaging method.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FX completed the conceptualization, investigation, methodology, data curation, formal analysis, writing-original draft, writing-review & editing in the manuscript. JL completed the Idea generation, methodology, data curation, writing-original draft, writing-review & editing in this manuscript. ZG completed the funding acquisition, data curation, formal analysis, writing-review & editing, supervision in this manuscript. JL completed the funding acquisition, formal analysis, writing-review & editing, supervision in this manuscript.