Anisotropic frequency-dependent characteristics of PP- and PS-waves in partially saturated double-porosity rocks


 Most frequency-dependent AVO inversions are currently based on an approximate equation derived using an isotropic medium. However, actual reservoirs usually show anisotropy, such as shale reservoirs, tight sandstone reservoirs and fractured reservoirs. We propose a joint frequency-dependent AVO (JFAVO) inversion in an anisotropic medium based on a periodic layered double-porosity medium. This JFAVO will allow us to quantitatively study the influence of fluids on the dispersion of PP- and PS-wave velocities and anisotropic parameters. First, we used a double-porosity medium to analyse the frequency-dependent characteristics of velocities and anisotropy parameters. We found that the anisotropic parameters show obvious dispersions, similar to those of velocities. Then, we derived the JFAVO inversion based on Rüger's equation to extract the dispersion of velocities and anisotropic parameters. Finally, we analysed the stability and applicability of the inversion algorithm, and used three sets of models to analyse the sensitivity of dispersion properties to fluids. The numerical analysis results show that PP-wave velocity dispersion and anisotropic parameter δ dispersion are sensitive to fluids, whereas, the velocity dispersion of the PS-wave is not. When saturation exceeds 80%, the velocity dispersion and anisotropic parameter dispersion properties are not sensitive to fluids.


Introduction
When seismic waves propagate in heterogeneous porous media, the pressure differences between the pores lead to the relative displacement between the fluid and solid phases. This phenomenon is called wave-induced fluid flow (WIFF), an important cause of seismic-wave dispersion. In seismic exploration, mesoscopic WIFF (much higher than pore scale and smaller than seismic wavelength) is the theoretical basis for quantitative predictions of reservoir fluids.
Seismic-wave dispersion is widely used in the fluid detection of layered porous reservoirs. White et al. (1975) put forward a layered porous medium filled with gas and brine. Norris (1993) used the homogenization method to discuss the influence of interlayer flow on P-wave dispersion.
On the basis of the propagation matrix theory, Brajanovski et al. (2005) calculated the P-wave modulus with a vertical incident in a double-porosity medium. Using a partially saturated model, Carcione & Picotti (2006) analysed the influence of the physical parameters of reservoirs on P-wave dispersion. They concluded that P-wave dispersion is a good fluid indicator. Xu et al. (2011) used an asymptotic expression to calculate the frequency-dependent reflection coefficient of fluid-saturated porous media and applied it in fluid identification in a gas field. However, these models are all based on the assumption of a vertical-incident P-wave and ignore S-wave velocity dispersion. Krzikalla & Müller (2011) derived the stiffness tensor for periodically layered porous media and analysed the P-and SV-wave velocity dispersion. Using this method, Kong et al. (2016Kong et al. ( , 2017 and Carcione et al. (2013) obtained five complex parameters for a transversely isotropic (TI) medium. Many scholars use the White theory to study the saturation of porous media (Lambert et al. 2006;Ren et al. 2009;Liu et al. 2011;Pang et al. 2016). Ba et al. (2011Ba et al. ( , 2016Ba et al. ( , 2017 proposed a doubleporosity model to study mesoscale attenuation in a poroelastic medium. Zhao et al. (2018) deduced the wave equation of a layered double-porosity medium model, and analysed the influence of fluid viscosity and permeability on attenuation and dispersion. Using the double-porosity medium, Yin et al. (2018) studied the frequency-dependent characteristics of the reflection and transmission waves of inhomogeneous media. They discussed the impact of rock parameters on the dispersion curves of reflection and transmission coefficients. Zhang et al. (2019) analysed the local squirt flow and global flow of porous media and studied the effects of different crack parameters on P-wave dispersion and attenuation. Barak et al. (2018) proposed a three-layered mathematical model, in which the middle layer was a double-porosity medium. Based on the double-porosity dual-permeability model, Kumar et al. (2018Kumar et al. ( , 2019 derived the expression of the P-wave reflection coefficient with angles and analysed the effects of frequency, fluid viscosity and double-porosity structure on the reflection coefficient. With the development of viscoelastic anisotropic media theory, the dispersive attribute gradually became an important reservoir fluid-detection parameter. Frequencydependent amplitude variation with offset (FAVO) is a useful method for fluid detection using seismic data, which has been used in oil and gas exploration for many years and achieved good results. Chapman et al. (2006) used the mesoscopic fracture model to analyse the attenuation and velocity dispersion of a reservoir filled with different fluids. They considered that velocity dispersion characteristics could identify reservoir fluid types, providing a theoretical basis for fluid detection using dispersive attributes. Wilson et al. (2009) proposed a FAVO inversion algorithm for the first time using the time-frequency analysis method and the Zoeppritz ap-proximate expression to extract velocity dispersion from seismic data. Hao et al. (2013) extended Shuey's approximate equation into the frequency domain to obtain a linear expression of frequency-dependent reflection coefficient, and then extracted the velocity dispersion of PP-wave using Bayesian theory. Zhao et al. (2014) comprehensively discussed the development process and trend of the FAVO technology, and considered that it could be widely used in quantitative reservoir prediction. Huang et al. (2017), Luo et al. (2018) and Wu et al. (2014) improved spectral decomposition to obtain a more accurate FAVO inversion. However, these inversion methods are based on PP-waves without considering anisotropy.
For the first time, we attempt to analyse the frequency variation characteristics of anisotropic parameters in doubleporosity media and deduce the joint frequency-dependent AVO ( JFAVO) inversion equation to extract the dispersion of PP-and PS-waves and anisotropic parameters. First, we calculate the elastic stiffness matrix of a partially saturated double-porosity model using the improved double-porosity medium model (Kong et al. 2013) and the scaling theory of elastic parameters (Krzikalla & Müller 2011). Then, we derive a JFAVO inversion method using Rüger's equation (Rüger 1997) to extract velocity dispersion and anisotropic parameter dispersion. Finally, we give three sets of numerical examples to analyse the effects of fluids on the four dispersions.

Layered double-porosity model
As we all know, sedimentary rocks occupy a large proportion of crustal rocks, and their most prominent feature is their layering characteristics. Porosity and permeability of sedimentary rocks usually change in the vertical direction. Hence, a layered double-porosity medium model can be used to analyse the gas saturation of real reservoirs quantitatively. Sandstone thin interbedded reservoir and shale gas reservoir can be simplified as a simple three-layer model. The middle layer represents the reservoir, as shown in figure 1a. Figure 1b is the schematic diagram of the layered double-porosity model, where layer 1 and layer 2 have different porosity and permeability.
When seismic waves propagate in layered double-porosity rocks, mesoscale interlayer flow is the main factor causing seismic wave attenuation. There are two steps in the construction of the model. First, we can calculate the modulus of P-wave for a layered double-porosity medium under the vertical incidence of the P-wave using the equation proposed by Brajanovski et al. (2005). Then, we can use the upscaling theory proposed by Krzikalla & Müller (2011) to obtain the elastic parameters of the double-porosity model. The gas saturation simulation is realized by changing the thickness of the gas-bearing layer. The effect of gas saturation on the dispersion of velocities and anisotropy parameters can be analysed on the basis of this model. Brajanovski et al. (2005) gave an expression for the P-wave modulus of a double-porosity medium. The model consisted of two groups of porous media with periodic distributions, and each period was composed of two isotropic layers (layer 1 and 2). Kong et al. (2013) proposed an expression of the P-wave modulus when the layered double-porosity medium was saturated with two kinds of fluid. The P-wave modulus consists of three terms

Stiffness tensor of double-porosity media
) . (2) Subscripts 1 and 2 represents layer 1 and 2, respectively. Ignoring the subscripts, the skeleton parameters include the bulk modulus K, shear modulus , permeability and porosity . The fluid parameters include bulk modulus K f , viscosity and density f . The Gassmann equation can be used to calculate the modulus of a P-wave.
where j represents Biot's coefficient.
Omitting the subscripts, Land M are the P-wave modulus and pore modulus of dry rocks, and Kand are the bulk and shear modulus of dry rocks obtained by the equation proposed by Krief et al. (1990): where represents porosity. Krzikalla & Müller (2011) used the elastic parameter upscaling theory to obtain the stiffness matrix C ij (i, j = 1,2, … ,6) of a vertically transversely isotropic (VTI) model. C ij depends on the stiffness matrixes of the rock in unrelaxed and relaxed stiffness state C u ij and C r ij , respectively. C u ij and C r ij can be calculated using the Backus model (Krzikalla & Müller 2011).
where R( ) is the scaling factor that can be calculated using the elastic parameters under vertical incidence of a P-wave, Once the elastic stiffness matrix of the model is obtained, we calculate the phase velocities using the equations proposed by Carcione et al. (2013).
The velocity and attenuation can be calculated using the following equations (Wang 2008).
where Re(M) and Im(M) represent the real and imaginary parts of the elastic modulus M, where M is elastic modulus, M = v 2 .

Joint frequency-dependent AVO inversion
The reflection coefficients on two elastic VTI surfaces can be obtained using Rüger's approximate expression assuming weak anisotropy and impedance (Rüger 1997). We arrange the Rüger approximate expression and use the following equations to obtain the PP-wave reflectivity and PS-wave reflectivity, where and are the Thomsen (1986) parameters.
Neglecting the effect of frequency on the density of media, Rüger's approximate equation can be rewritten in the frequency domain as follows: Equations (16) and (17) are expanded using the Taylor series at the reference frequency f 0 , and the first-order terms are retained, which have the following expressions: where I a and I b are the velocity dispersions of PP-and PSwaves at frequency f 0 , and I c and I d are the dispersions of anisotropic parameters at frequency f 0 . Velocity dispersion is the derivative of the velocity reflection coefficient to a frequency and the dispersion of anisotropic parameters represent the derivatives of the difference of anisotropic 358 parameters to a frequency.
The matrix expression of the PP-and PS-waves is The equations of joint frequency variant AVO inversion can be expressed as Omitting the lower subscripts PP and PS, dR is the reflection coefficient difference between the reference and the inversion frequencies, Sis the forward coefficient matrix and m is the dispersion attribute to be solved. If the frequency involved in the inversion is m and incidence angle is n, then the joint inversion equations are composed of 2(m − 1)n equations. This is an overdetermined problem, and we can use the regularization method to solve it.
For JFAVO inversion, we can construct the following objective function: where is the weight coefficient that indicates the proportion of the PP-wave data in the joint inversion. We seek the partial derivative of the objective function If the partial derivative is zero, the model solution of the joint frequency variant AVO inversion can be obtained.

Numerical simulation and results
The periodically layered porous medium is filled with gas and brine. We assume that the mineral grains of the rock are quartz, with parameters K g = 38GPa, g = 44GPa and g = 2.65kg∕m 3 . The periodic thickness of the model is 1 m, and the porosities of layers 1 and 2 are 30 and 20%, respectively. The permeability can be calculated by the Kozeny-Carman relation = B 3 D 2 (1− ) 2 (Mavko et al. 2009), where D is the grain diameter and B = 0.003. The bulk modulus K f , viscosity and density f of fluid are given in Table 1.

Influence of fluid distribution characteristics on double-porosity media
Four models are designed to compare and analyse the velocity dispersion of P-and SV-waves in different distribution characteristics of fluids. The periodic thickness is 1 m, and the thickness of layers 1 and 2 is 0.5 m. Table 2 shows the fluid distribution characteristics of models 1-4. Models 1 and 4 are completely filled with water and gas, and the corresponding gas saturations are 0 and 100%, respectively. For model 2, layers 1 and 2 are filled with water and gas, respectively, and the gas saturation is 40%. In model 3, layer 1 is filled with gas and layer 2 is saturated with water. The gas saturation is 60%. Using these equations, we can obtain the velocities, anisotropy parameters and attenuation coefficients of the four models. Figure 2 shows the velocities of P-and SV-waves and their corresponding attenuations when the incident slowness angle is 30°. In the figure, the blue line corresponds to model 1, the red line to model 2, the pink line to model 3 and the green line to model 4. From figure 2, we can conclude that the velocities of model 4 do not obviously change with frequency. Figure 2 parts a and b show that the P-wave velocities of models 1, 2 and 3 increase with increasing frequency, and the P-wave velocity of model 2 increases most obviously. The attenuation coefficients of P-wave first increase and then decrease with the increasing frequency. The attenuation of model 2 is the highest, and its corresponding peak frequency is about 100 Hz.  coefficient of SV-wave first increases and then decreases with increasing frequency, and the peak value of the attenuation coefficient of model 2 is the largest. Within the seismic frequency band (0-100 Hz), the attenuation coefficients of shear wave from large to small follow the order of models 1, 2, 3 and 4. Figure 3 presents the frequency-dependent characteristics of the anisotropy parameters and . From figure 3a, it is clear that the anisotropy parameter of models 1 and 2 decreases as frequency increase. In contrast, the anisotropy parameter of model 3 increases with increasing frequency and that of model 4 is constant at 0.23. Figure 3b shows the variation characteristics of the anisotropy parameter with frequency. Except for model 4, the anisotropy parameter increases with increasing frequency. The value of the anisotropy parameter of model 4 is small and remains near zero.
After obtaining the elastic stiffness matrix of the doubleporosity medium, we use the Zoeppritz equation rederived by Qin et al. (2018) to calculate the reflection coefficients of the surface between an isotropic medium and a doubleporosity medium. Because the PP-wave dispersions of model 1 and 4 are not evident, for simplicity, we only show the reflection coefficients of PP-and PS-waves for models 2 and 3. Figure 4 parts a and b show the relationship between the reflection coefficients and the incident angle of model 2. Different colours of the curves in the figure represent different frequencies. The absolute values of the PP-wave reflection coefficients decrease with the increase of frequency, and decrease with increasing slowness angle. Figure 4 parts c and d show the reflection coefficients of model 3. Comparing the reflection coefficients of models 2 and 3, we found that the divergence of the reflection coefficients for model 2 is significantly higher than that of model 3, and the divergence of PPwave reflection coefficients is obviously higher than that of the PS-wave.
JFAVO inversion is a linear problem and its expression in equation (24) can be simplified as where dis frequency-dependent reflectivity, Gis the forward operator and m the unknown parameters of dispersion. Gis completely dependent on the model parameters and velocity field. The stability of the JFAVO inversion algorithm can be analysed by the condition number of G (Wang 1999;Downton 2005). Figure 5 shows the condition numbers as functions of the maximum angle used in JFAVO for model 2. The curves with different colours correspond to different weight coefficients, representing the proportion of PPwave involved in the inversion. When the maximum angle is greater than 20°, the condition numbers of different weight coefficients are less than 100 and the inversion result is stable. Therefore, this method can be applied to seismic field data.
Using this inversion method, we performed the JFAVO inversion on the reflection coefficients PP-and PS-waves, and obtained the velocity dispersion of PP-and PS-waves and dispersion of anisotropic parameters. We selected the reflection coefficients of 10, 30, 50, 70 and 90 Hz as the input parameters, with 50 Hz as the reference frequency. We invert the reflection coefficient of the theoretical model data and the reflection coefficient with 5% random noise to analyse the stability of the inversion algorithm. Figure 6 shows the results. In figure 6a, the blue and red curves represent the velocity dispersion of PP-and PS-waves, and the pink and green curves are the dispersion of the anisotropic parameters and , respectively. The abscissas 1-4 represent models 1-4, respectively. It is evident from figure 6 that the dispersion of PP-wave velocity and anisotropic parameters of model 2 are significantly higher than those of the other three models, and the dispersions of PS-wave velocity of the four models are close to zero. Figure 6b shows the inversion results of the theoretical model with 5% random noise, which is consistent with the inversion results of the theoretical models.

Double-porosity model based on Brie law
By comparing models 1 and 4, we found that the velocity dispersions of PP-and PS-waves for the double-porosity model are different when it is saturated with gas and water. For analysing the effect of fluid properties on the dispersions of the model, we assume that the double-porosity sandstone reservoir is filled with uniformly mixed fluid. The Wood equation is mainly applicable to isotropic porous media, so we use Brie's law for fluid mixing. Brie et al. (1995) gave an equation for obtaining a gas-water mixture's equivalent bulk modulus.
where K f , K b and K g express the bulk modulus of the mixed fluid, brine and gas, respectively; S b expresses water saturation and e is an empirical parameter. Fang et al. (2019) analysed the influence of pore fluid saturation on PP-wave velocity through laboratory measurements and confirmed the applicability of Brie's equation. According to equation (29), the bulk modulus of the mixed fluid increases with increasing water saturation. By introducing bulk modulus of the mixed fluid into the calculation equation of the double-porosity medium model, the characteristics of P-and SV-wave velocity and anisotropy parameters varying with frequency under different gas saturations can be obtained. Figure 7 parts a and b show the velocities and attenuation of a P-wave. Figure 7 parts c and d show the velocities and attenuation of an SV-wave. The velocities of P-and SV-waves increase with increasing frequency, and the maximum attenuation decreases with the increasing gas saturation. The frequency corresponding to the peak attenuation decreases with increasing gas saturation.  different gas saturations. The anisotropy parameter decreases with increasing frequency and increases with increasing gas saturation. As the frequency increases, the attenuation first increases and then decreases. As the gas saturation increases, the peak attenuation gradually decreases. The attenuation of the anisotropy parameter drops to about zero when the gas saturation is 100%. Because the frequency characteristics of the anisotropy parameter are similar to those of , the curves about are not shown for simplicity. The dispersion and attenuation caused by mesoscopic interlayer fluid flow increases with the increase of bulk modulus.
Through JFAVO inversion, velocity dispersions of PPand PS-waves and dispersions of the anisotropic parameters and are calculated, which are functions of gas saturation shown in figure 9. The results of the JFAVO inversion show that the four dispersions decrease with increasing gas saturation, and the minimum value is obtained when the model is completely saturated with gas. The velocity dispersion of PPwave is similar to that of anisotropic parameter , whereas, the velocity dispersion of PS-wave is the smallest, which is about zero. The anisotropic parameter dispersion is higher 362 Figure 6. Velocity dispersion of PP-(blue) and PS-waves (red), dispersion (pink) and dispersion (green) of the theoretical model data (a) and data with noise (b) for the four models.  than those of others, which is approximately twice as much as the dispersion of PP-wave. When the gas saturation is less than 50%, the four kinds of dispersion are approximately linear relative to gas saturation. The four kinds of dispersion decrease slowly with increasing gas saturation when the saturation is greater than 50%.

Saturated gas and brine in double-porosity sandstone reservoirs
We assumed that layers 1 and 2 are filled with different fluids to study the double-porosity medium saturated with two immiscible fluids. With a fixed periodic thickness H, different saturations can be simulated by adjusting the thickness of layer 1 and layer 2. In this part, we designed two models. Model A is saturated with water in a high-porosity formation, that is, layer 1 is filled with water. Model B is saturated with gas in high-porosity formation, that is, layer 1 is filled with gas.
Using these equations, we can calculate the frequencydependent velocities and attenuation. When the incident angle is 30°, the velocities and attenuation of model A are plotted as functions of the frequency with different gas saturations in figure 10. The P-wave velocities and attenuation are shown in figure 10 parts a and b, and those of an SV-wave are shown in figure 10 parts c and d, respectively. When the gas saturation is less than 80%, the P-wave velocities have an obvious dispersion, whereas the velocity of the SV-wave has no obvious dispersion. The velocities of P-and SV-waves decrease with increasing gas saturation. When gas saturation increases by 10% from 0.1 to 1, the peak value of the P-wave attenuation coefficient decreases with increasing gas saturation. In contrast, the peak attenuation of SV-wave increases first and then decreases with increasing gas saturation. The maximum attenuation occurs at 60%. The attenuation coefficient of the SV-wave is obviously smaller than that of the P-wave. Figure 10 parts e and f show the frequencydependent anisotropy parameter and the attenuation coefficient with different gas saturations. The anisotropy parameter decreases with increasing frequency, and first increases and then decreases with increasing gas saturation. The anisotropy parameter is the largest with 40% gas saturation. With the increase of frequency, the attenuation of first increases and then decreases. Except for when the gas saturations are 0 and 100%, the peak attenuation of is almost equal and the corresponding frequency increases with increasing gas saturation. When the gas saturations are 0 and 100%, the model becomes a homogeneous isotropic medium saturated with water and gas, respectively. In this case, the anisotropy parameter and attenuation coefficient of the model are both zero. 364 Figure 10. Velocities of (a, b) P-and (c, d) SV-waves, (e, f) anisotropy parameter and attenuation coefficients vary with frequency of model A.
When the incident angle is 30°, the velocities and attenuation coefficients of model B are drawn as functions of the frequency with different gas saturations in figure 11. Figure 11 parts a and b show the velocities and attenuation of a P-wave. Figure 11 parts c and d show the velocities and attenuation of an SV-wave. The variation characteristics of P-and SVwave velocities and attenuation coefficients with frequency and gas saturation of model B are similar to those of model A. At the same gas saturation, the peak values of the attenuation coefficients of PP-and SV-wave of model B are smaller than those of model A. The corresponding frequency of the peak values is also smaller than that of model A. Figure 11 parts e and f show the anisotropy parameter and the attenuation coefficient as frequency functions with different gas 365 Figure 11. The velocities of (a, b) P-and (c, d) SV-waves, (e, f) anisotropy parameter and attenuation coefficients vary with frequency of model B.
saturations. The anisotropy parameter increases as the frequency increases, which is contrary to that of model A. The anisotropy parameter is the largest, with 60% gas saturation. When the double-porosity medium is filled with gas and water, the peak values of attenuation coefficients are almost the same and positive. Figure 12 shows the JFAVO inversion results of the two models with different gas saturations. Figure 12 parts a and b correspond to the inversion results of models A and B, respectively. From figure 12a, the dispersion of PP-wave velocity is greater than that of the anisotropic parameters and the dispersion of PS-wave is around zero. When the gas saturation is 1-5%, the four dispersions increase rapidly and then decrease as the gas saturation increases. When the gas saturation is greater than 60%, the velocity dispersions of the PPwave and the dispersion of tend to zero. In contrast, the 366 Figure 12. Velocity dispersions of PP-(blue) and PS-waves (red) and dispersions of the anisotropy parameters (pink) and (green) of models A (a) and B (b).
dispersions of PP-wave velocity and anisotropic parameter tend to zero with gas saturations greater than 40%. Figure 12b shows that the dispersions of PP-wave velocity and the anisotropy parameter first increase and then decrease with increasing gas saturation. In contrast, the velocity dispersion of the PS-wave is around zero. The dispersion of anisotropy parameter presents a bimodal shape with gas saturation. When gas saturations are greater than 80%, the values of the four dispersions are around zero. When gas saturations are lower than 7%, the velocity dispersion of PP-wave increases rapidly with increasing gas saturation and reaches the maximum. When gas saturation is between 7 and 80%, the velocity dispersion of PP-wave decreases with increasing gas saturation. With gas saturations greater than 80%, the dispersion of PP-wave tends to zero. Overall, the velocity dispersion of PS-wave is obviously smaller than other dispersions, which is about zero. When gas saturations are lower than 5%, the dispersion of anisotropy parameter increases sharply. When the gas saturation is between 5 and 40%, the dispersion of anisotropy parameter increases slowly with increasing gas saturation. When the gas saturation increases from 4% to 80%, the dispersion of anisotropy parameter decreases rapidly. When the model is saturated with gas and water, the velocity dispersion of PP-wave and attenuation caused by mesoscale interlayer fluid flow first increase, then decrease with gas saturation and reach the maximum value when gas saturation is small.

Conclusions
We analysed the frequency variation characteristics of anisotropic parameters in a double-porosity medium and derived the JFAVO inversion equation to analyse the impact of fluids on the velocity dispersions of PP-and PS-waves and the dispersions of anisotropic parameters. The numerical analysis results show that the influence of fluids on the velocity dispersion of a PS-wave is small. When the double-porosity medium is saturated with a single fluid or a uniformly mixed fluid, the velocity dispersion of a PP-wave and the dispersion of anisotropic parameters decrease with decreasing gas saturation. However, when the double-porosity medium is saturated with two-phase immiscible fluids, the velocity dispersion of a PP-wave and dispersion of anisotropic parameter first increase and then decrease as gas saturation increases. The dispersion of is most sensitive to saturation when the high-porosity layer is saturated with gas. The dispersion of a PP-wave is most sensitive to saturation when the high-porosity layer is saturated with water in the model. When the gas saturation is greater than 80%, the dispersion properties of a double-porosity medium do not change obviously with gas saturation.