Numerical Investigation of the Effect of Heterogeneous Pore Structures on Elastic Properties of Tight Gas Sandstones

Strong heterogeneity of pore microstructures leads to complicated velocity-porosity relationships in tight sandstone that cannot be well explained by conventional empirical formulas. To better understand the effect of complex pore structures on elastic properties of tight gas sandstone, we compared three rock physics models. In the first model, we used a single aspect ratio value to quantify varied pore geometry in the tight sands. In the second model, complex pore space was equivalent to the combination of high-aspect-ratio round pores (stiff pores) and low-aspect-ratio compliant microcracks (soft pores). In the third multiple pore-aspect-ratio model, pore spaces are represented using a set of pores with varied values of aspect ratio following statistical normal distribution. Modeling results showed that complex velocity-porosity relationships could be interpreted by the variations in pore aspect ratio in the first model, by the fraction of soft pores in the second model, and by the mean value and variance in the third statistical model. For a given mean value in the third model, higher variance of the multiple pore-aspect-ratio indicated stronger heterogeneity of pore spaces. Further studies on rock physical inversion showed that, compared with the first single pore-aspect-ratio model, the second dual-pore model gave better prediction in shear wave velocity by regarding the soft pore fraction as a fitting parameter. This finding revealed that the dual-pore model could be a more realistic representation of tight sandstone. The third statistical model showed comparable precision in the prediction of shear wave velocity compared with the dual-pore model; however, uncertainty existed for simultaneously determining mean value and variance of pore aspect ratio. On the basis of the dual-pore model, we evaluated the elastic modulus of dry frames of the tight sandstone using logging data in a borehole. Compared with empirical formulas, such as the Krief methods, the method in this paper provided a more rigorous way to determine elastic properties of dry frames for the tight sandstone. Comparisons of rock physical modeling methods offer a better understanding of the microstructures controlling the elastic behaviors of tight gas sandstone.

Strong heterogeneity of pore microstructures leads to complicated velocity-porosity relationships in tight sandstone that cannot be well explained by conventional empirical formulas. To better understand the effect of complex pore structures on elastic properties of tight gas sandstone, we compared three rock physics models. In the first model, we used a single aspect ratio value to quantify varied pore geometry in the tight sands. In the second model, complex pore space was equivalent to the combination of high-aspect-ratio round pores (stiff pores) and low-aspect-ratio compliant microcracks (soft pores). In the third multiple pore-aspect-ratio model, pore spaces are represented using a set of pores with varied values of aspect ratio following statistical normal distribution. Modeling results showed that complex velocity-porosity relationships could be interpreted by the variations in pore aspect ratio in the first model, by the fraction of soft pores in the second model, and by the mean value and variance in the third statistical model. For a given mean value in the third model, higher variance of the multiple pore-aspect-ratio indicated stronger heterogeneity of pore spaces. Further studies on rock physical inversion showed that, compared with the first single poreaspect-ratio model, the second dual-pore model gave better prediction in shear wave velocity by regarding the soft pore fraction as a fitting parameter. This finding revealed that the dual-pore model could be a more realistic representation of tight sandstone. The third statistical model showed comparable precision in the prediction of shear wave velocity compared with the dual-pore model; however, uncertainty existed for simultaneously determining mean value and variance of pore aspect ratio. On the basis of the dual-pore model, we evaluated the elastic modulus of dry frames of the tight sandstone using logging data in a borehole. Compared with empirical formulas, such as the Krief methods, the method in this paper provided a more rigorous way to determine elastic properties of dry frames for the tight sandstone. Comparisons of rock physical modeling methods offer a better understanding of the microstructures controlling the elastic behaviors of tight gas sandstone.

INTRODUCTION
Tight sandstone of low porosity and low permeability usually represent complex pore structures, which significantly impacts elastic properties of the tight sandstone. Complicated elastic behaviors of the tight sandstone with heterogeneous pore structures cannot be properly described by traditional rock physical methods that rely on empirical relations. Thus, it is necessary to investigate elastic properties related to heterogeneity using rock physical methods that incorporate pore structures.
The effect of microcracks on the rock elastic properties has been studied by Stuart and David (1977) in their theoretical rock physics models. After that, numerous works have been dedicated to the connection between elastic properties and microstructures in the tight sandstone reservoirs. According to study by Tutuncu et al. (1994), the existence of microcracks with low aspect ratio leads to the phenomenon that the P-and S-wave velocities Vp and Vs are sensitive to pressure changes. The presence of low-aspect-ratio pores or microcracks can explain the complex relationship between velocity and porosity in tight sandstone gas reservoirs (Smith et al., 2009;Vernik and Kachanov, 2010). Ruiz and Cheng (2010) introduced two sets of pores with different aspect ratios in the rock physical modeling, and applied the model to the prediction of elastic modulus of a tight sandstone. Yin et al. (2017) conducted experimental measurements and rock physical modeling to investigate the effect of pressure and fluid on elastic moduli of the tight sandstone, and found that microcracks in tight sandstone dominate the changes in elastic properties of the rock. Yan et al. (2011) employed the Mori-Tanaka model (Mori and Tanaka, 1973) to analyze experimental measurements on the dry frame of sandstone of low porosity. Deng et al. (2015) explored the influence of pore structures on velocity dispersion of tight sandstone, and found that the existence of microcracks is related to the heterogeneity of the rock skeleton and the fluid distribution. Ba et al. (2017) established a double dual-pore model for the anelasticity of P-wave propagation to analyze the effect of dissipation mechanisms caused by the heterogeneity of fabric and fluid distribution. Wang (2017) proposed a rock physics model for the prediction of gas-bearing sandstone properties based on the measurements and analysis on data for tight sandstone samples. Wang et al. (2020) incorporated the Bayes discriminant and deterministic rock physics modeling method to implement joint probabilistic fluid discrimination of tight sandstone reservoirs.
In this study, we built three rock physics models for the description of elastic properties associated with heterogeneous pore structures in tight gas sandstone. Using the three built models, we then predicted shear wave velocity and inverted corresponding parameters of pore structure to evaluate and compare the applicability of the models. Finally, we used the built model to estimate elastic modulus of the rock skeleton of tight gas sandstone.    Figure 1 shows the logging data from the tight gas sandstone reservoir in the Ordos Basin, China. The strata consist of the interbed of sandstone and mudstone. The tight sandstone corresponds relatively local low value of GR. Compared with the mudstone, the sandstone has relatively higher P-wave and S-wave velocities around 4.5 km/s and 2.5 km/s, and a relatively lower density of approximately 2.5 g cm −3 . Porosity generally is below 15% for the sandstone. As illustrated by the logging data of water saturation, the sandstone has a gas saturation usually higher than 40%. Data of the mineralogical volume show variations in the fraction of quartz and clay in the formation.

WELL DATA OF TIGHT GAS SANDSTONE RESERVOIR
Using the logging data in Figures 1, 2 show the variations in P-wave velocity Vp and S-wave velocity Vs versus porosity φ for sandstone and mudstone in the strata. Although velocities generally decrease with increasing porosity, the data points represent certain degree of scatting distributions, which indicates that porosity is not the only factor that affected the elastic properties of the tight sandstone. According to the discussion given by Smith et al. (2009) and Ruiz and Cheng (2010), both the volume of pore space (porosity) and the distribution or geometry of the pore (such as aspect ratio) have an impact on the elastic modulus of the tight sandstone. Thus, it is necessary to establish rock physical models that consider complex pore structures to quantitatively predict elastic behaviors for tight sandstone reservoirs.

ROCK PHYSICS MODEL Theory
As illustrated by the schematic diagrams in Figure 3, we compare three methods for modeling heterogeneous pore structures of tight sandstone. Figure 3A shows the single pore-aspect-ratio model, where all pores have the same geometry represented by the pore aspect ratio. Figure 3B shows the dual-pore model, where the total porosity φ Total is equivalent to the combination of two parts, including the stiff pores with a porosity φ Stiff and soft pores with a porosity φ Soft , that is, The stiff pores are set to be spherical and have an aspect ratio α Stiff = 1, whereas the soft pores are assumed to have an aspect ratio α Soft = 0.01. We use the fraction of soft pores in the total pore space, that is, f Soft = φ Soft / φ Total , to represent heterogeneous distribution of the pore space in the dual-pore model. Figure 3C illustrates that the value of the pore aspect ratio is assumed to have the normal distribution. In this model, heterogeneous distributions of pore space are controlled by the mean value α M of the pore aspect ratio and corresponding variance α V . Accordingly, Figure 4 shows two cases of normal distribution of pore aspect ratio, where α M and α V are 0.4 and 0.01 for the first scenario, and 0.6 and 0.0001 for the second. A higher value of α M means that the pores tend to be round and a higher α V implies that the pore geometry has a wide range of shapes.
The rock physical methods for the three cases in Figure 3 are illustrated in Figure 5. In all the three models, we use the average of the Hashin-Shtrikman bounds (HSB) (Hashin and Shtrikman, 1963) to obtain elastic properties of solid matrix, and then apply the self-consistent approximation (SCA) method (Berryman, 1995;Mavko et al., 2009) to add pores with different geometries and distributions as discussed above. Elastic properties of minerals and fluids used in the rock physical model are given in Table 1. FIGURE 6 | (A) P-wave velocity Vp and (B) S-wave velocity Vs varying with pore aspect ratio α and porosity φ for the single pore-aspect-ratio model of the tight sandstone. (C) P-wave velocity Vp and (D) S-wave velocity Vs varying with porosity φ for a series of pore aspect ratio α for the single pore-aspect-ratio model with logging data from the target zone overlapped.
Frontiers in Earth Science | www.frontiersin.org For the construction of the rock physical models in Figure 5, the SCA theory for multi-phase media is as follows: where j is the j-th material, f j is its volume fraction, K j and µ j represent the bulk modulus and shear modulus of the j-th component, K * SC and µ * SC represent the calculated equivalent bulk modulus and shear modulus, and P * i and Q * i are geometric factors of the inclusion added to the background medium.
The bulk module of the fluid mixture consisting of gas and water is calculated by the Wood (1955) theory, as follows: where K w and K g are bulk modulus of water and gas, respectively, S w and S g are water and gas saturation, respectively.
The density of the fluid mixture is as follows: where ρ w and ρ g are the densities of water and gas, respectively. Figures 6A,B show the Vp and Vs varying with pore aspect ratio α and porosity φ for the single pore-aspect-ratio model of the tight sandstone. It can be observed that Vp and Vs decrease rapidly with increasing porosity, especially for lower-valued pore aspect ratio α, because the presence of elongated pores can dramatically reduce the elastic modulus of rocks. Accordingly, Figures 6C,D illustrate calculated Vp and Vs varying with porosity φ for a set of pore aspect ratio α that has a value ranging from 0.05 to 0.2. Corresponding logging data from the target zone are overlapped. The relationship between Vs and φ is more scattered compared with the relationship between Vp and φ.

Rock Physical Modeling
In Figures 7A,B, the presence of small amount of soft pores can dramatically reduce velocities of the tight sandstone. Accordingly, Figures 7C,D indicates that the parameter f Soft ranges from 0.05 to 0.2 for the Vp data and from 0.05 to 0.4 for the Vs data for the interpretation of the velocityporosity relationships.
Similarly, Figures 8A,B show velocities that change for varied mean value α M and variance α V of pore aspect ratio in the multiple pore-aspect-ratio model. A higher value of α M corresponds to the presence of more spherical pores and thus results in higher modulus and velocities. The value of α V denotes the spread range of the pore aspect ratio, and  a higher α V means more heterogeneity of pore geometry distribution. For a specific value of α M , the increase in α V reduces velocities.
In Figures 8C,D, we assume α M = 0.75 in the tight sandstone. We find that α V should vary from 0.01 and 0.06 to interpret the velocity-porosity relationships of the overlapped logging data. Figures 6-8 show rock physical modeling using the constructed models in Figure 5, where the procedures for calculation are described in details. The purpose of the rock physical modeling is to illustrate how the porosity and pore geometry parameters influence the Vp and Vs. As shown in Figures 6-8, it means that the model can interpret the data if the overlapped data fall into the ranges covered by the calculated curves. Then, it is feasible to conduct rock physical inversion to estimate pore geometry parameter as shown in Figure 9.

Vs PREDICTION USING THE BUILT ROCK PHYSICS MODELS
On the basis of the constructed rock physics models shown in Figure 5, we predict S-wave velocity Vs, when the information of porosity, mineral compositions, and fluid saturation is known as illustrated in Figure 1.
According to the method by Guo and Li (2015), we find the best fit parameters, including pore aspect ratio α, the proportion of soft pores f soft , and the variance of pore aspect ratio α V in the three models by adjusting the values of them to minimize the absolute difference between the P-wave velocities Vp measured in the well and those calculated by the rock physics models in Figure 5. The S-wave velocities Vs calculated by the estimated α, f soft and α V are regarded as the predicted values. Figures 9A-C show the prediction results from the three rock physics models, where all of the prediction results provide reasonable estimates of Vs. The correlation coefficients between the predicted and measured Vs are 0.9077, 0.9375, and 0.9277, respectively. As illustrated in Figure 9B, the inversion based on the dual-pore model generates the best result for the Vs prediction. Theoretically, the multiple pore-aspect-ratio model in Figure 5C may work better, as there always exists a certain normal distribution of the multiple pore-aspect-ratio in the model to fit the true scenario. However, it may takes too many parameters to determine, which is unlikely to realize. Comparatively, the dual-pore model shows the advantage to obtain a relatively fine and stable result without solving too many parameters.

Calculation of Elastic Modulus for Dry Frame
The elastic modulus of rock skeleton that are determined by minerals, porosity, and pore structures are the key input parameters for various rock physical theories, including the Gassmann fluid substitution theory (Gassmann, 1951), the Biot theory (Biot, 1955) and the BISQ theory (Dvorkin and Nur, 1993;Dvorkin et al., 1994).
Although the elastic modulus of the rock skeleton can be obtained through laboratory measurements of specimens, it is difficult to estimate these properties for in situ underground formations. Empirical formulas usually are used to estimate elastic modulus of the dry frame.
Based on the experimental data from Raymer et al. (1980); Krief et al. (1990) proposed the nonlinear empirical relationships between the rock skeleton modulus and porosity, as shown by the followings: where K 0 and µ 0 are the bulk and shear modulus of the solid matrix, respectively, and ϕ is porosity. Empirical parameter m generally ranges from 0 to 10, which is related to porosity, pressure, temperature, clay content, and other lithology factors (Han et al., 1986). Because these empirical relationships between elastic modulus of rock skeleton and porosity were established based on data from specific study areas, it may limit the application of the methods to other field data which represent different lithology or microstructures.
In this work, based on the inverted results of pore structures using the dual-pore model shown in Figure 9B, we can calculate the elastic modulus of rock skeleton in a manner that represents more physical meaning. We conduct rock physical modeling to calculate the elastic modulus of rock skeleton by using the inverted values of pore structure parameter f soft , and by setting the bulk and shear modulus of the pore fillings to 0. The calculation results are shown in Figure 10.
Data points in Figure 11 show the elastic modulus of skeleton predicted by the method presented discussed above. Curves in Figure 11 are the elastic modulus estimated by the Krief empirical formulas in Eqs 5 and 6. Due to the presence of complex pore structure in tight sandstone, the relationship between the modulus of the skeleton and porosity are scattering and thus cannot be well described by adjusting the empirical parameter m in the Krief theory.
Comparatively, the methods presented in this study have considered the influence of minerals, porosity, and pore structures in the calculation of the elastic modulus of dry frame, which show more physical meaning compared to the Krief empirical relations.  data in Figure 1. In Figure 12A, the Poisson's ratio of tight sandstone is mainly between 0.2 and 0.3 when the porosity is less than 0.1. Porosity that is greater than 0.1 corresponds to higher gas saturation and relatively lower Poisson's ratio. Figure 12B is the crossplot of the Poisson's ratio versus porosity of the rock skeleton of the same sandstone, with calculated elastic modulus of the rock skeleton shown in Figure 10. In Figure 12B, it is found that the Poisson's ratio of the rock skeleton of the sandstone decreases with increasing porosity.
For further analysis, Figures 12C,D are the crossplot of the Poisson's ratio versus the micro-crack porosity φ soft (φ soft = φf soft ) of the fluid-saturated rock and the rock skeleton of the tight sandstone, respectively. Figure 12D shows that the micro-crack porosity has a linear relationship with the Poisson's ratio of the rock skeleton. The decrease in the Poisson's ratio of the dry frame corresponds to an increase in the fraction of microcracks. It reveals that the Poisson's ratio of the rock skeleton is sensitive to the heterogeneity of pore spaces.

DISCUSSION AND CONCLUSION
For the tight gas sandstone with low porosity and permeability, porosity is not the only factor that affects elastic behaviors of the sandstone. Heterogeneity of pore spaces should be considered in rock physical modeling to reasonably predict elastic modulus.
We established three rock physics model to describe the complex pore structure of tight gas sandstone reservoirs. In the single pore-aspect-ratio model, all pores are assumed to have the same pore aspect ratio. In the dual-pore model, the pore space is equivalent to a combination of stiff and soft pores. In the multiple pore-aspect-ratio model, the aspect ratio mean value and variance are used to represent statistical normal distribution of pores with varied aspect ratio.
Results for pore structure inversion and S-wave velocity prediction show that the dual-pore model works well and behaves as a more appropriate representation of the tight gas sandstone compared to the other model built in this study. The inverted parameter of soft pore fraction f Soft can explain the heterogeneity of pore structure of the tight sandstone. Increasing fraction of soft pores can dramatically decrease elastic modulus of the tight sandstone.
We applied the dual-pore rock physics model to estimate elastic modulus of rock skeleton. Compared to conventional empirical formulas, such as the Krief methods, the methods presented in this study considered the influence of minerals, porosity and pore structures in the calculation of the elastic modulus of dry frame of tight sandstone, which show more physical meaning and is more appropriate for the dry frame representing heterogeneous pore structures.
As to the applicability of the built rock physics models, it is not sure whether these models are valid for the tight sandstone from different areas. However, the rock physics models that are constructed using effective medium theories generally work well than those empirical models. The reason is that the former ones follow certain kinds of physical laws by incorporating parameters with specific physical meaning (for example, aspect ratio of pores). If the geological scenario agrees with the assumptions of the effective medium theories, the built model will work well and generate good results. The application results show that the models in this study are suitable for the tight sandstone of the research area. However, it may be that not all the tight sands are the same. So, the applicability of the built rock physics models should be carefully checked when applied to other fields.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
ZG contributed as the major author of the manuscript. XQ did a part of writing and coding works. YZ, CN, DW, and YL provided some interesting ideas. YZ provided the suggestions. CN and DW gave analysis on results. YL collected the data. All authors contributed to the article and approved the submitted version.