Estimation of the Pore Microstructure of Tight-Gas Sandstone Reservoirs with Seismic Data

Tight-sandstone reservoirs have a complex pore structure with microcracks and intergranular pores, which have a significant impact on the seismic properties. We have performed ultrasonic measurements at different confining pressures for 15 tight-gas sandstone samples of the Xujiahe formation in Western Sichuan Basin, and have available well-log and seismic data of this area. The aim of this work is to estimate the porosity and crack properties for variable pressure conditions. The EIAS (equivalent inclusion-average stress) model is adopted to compute the high- and low-frequency bulk and shear moduli as a function of crack aspect ratio and (soft) and (stiff) porosities. Then, we use the EIAS-Zener anelastic model to obtain the wave properties as a function of frequency, and compare results with those of the constant Q (Kjartansson) one for verification of the robustness of the approach. The corresponding P-wave impedance, density and phase velocity ratio (V P/V S) are computed in order to built 3D rock‐physics templates (RPTs) at the ultrasonic, well-log and seismic frequency bands. The methodology is applied to a survey line crossing two wells, which together with the laboratory experiments, provide calibration suitable data. The estimated stiff porosity and crack porosity and density are consistent with the available data and actual production records, indicating that 3D RPTs provide a useful interpretation tool in seismic exploration and prospect evaluation.


INTRODUCTION
Unconventional gas resources offer significant gas production growth potential. Tight-gas sandstones represent an important part of the unconventional production and abundant reserves are yet to be developed (Khlaifat et al., 2011). There are tight-gas sandstone reservoirs almost in all oil and gas bearing areas with significant gas reserves (Zhu et al., 2008). These reservoirs have sandstones with low porosity, permeability and gas saturation and abundant microcracks (Guan andNiu, 1995, Anjos et al., 2003;Cheng et al., 2019;Guo et al., 2018). An effective reservoir characterization is the key to their successful exploration (Storker et al., 2013;Liu et al., 2019). The reservoir contains a large number of natural fractures, which provide main seepage channels and effective storage space, and favor single-well productivity (Harmelen and Weijermars, 2018;Xiao et al., 2019). Therefore, an accurate identification of the pore properties is essential.
In this case, rock physics is required to relate the wave properties to the composition and microstructure of these rocks, which have been affected by complex diagenetic and sedimentary processes during their formation (Vernik and Kachanov, 2010;Zhang et al., 2019;Dvorkin et al., 2020). Walsh (1965) developed some of the first techniques to study the effects of cracks and confining pressure on the dry-rock moduli. Thomsen (1995) considered aligned circular cracks and analyzed the effects of the associated anisotropy on the elastic properties, and Smith et al. (2009) quantified how microcracks affect the wave velocities in tight-gas sandstone reservoirs. Tang (2011) proposed an elastic-wave model to describe the effects of pores and cracks by extending the Biot and BISQ theories. Double-porosity models that consider two phases with different compressibilities are relatively new and useful approaches to model wave propagation in heterogeneous media (Pride et al., 2004;Ba et al., 2016;Ba et al., 2017;Fu et al., 2018).
To model the relation between pore geometry and elasticity attributes, we used the EIAS (equivalent inclusion-average stress) model (Endres and Knight, 1997), which is consistent with Gassmann equation at low frequencies and with the Hashin-Shtrikman bounds when applied to two-phase systems. Moreover, we used the EIAS-Zener model that generalizes the EIAS model to the full frequency range, by incorporating the Zener (standard linear solid) model . The crack porosity, average crack aspect ratio and quality factors and phase velocities were obtained as a function of frequency by fitting experimental data at the unrelaxed and relaxed states.
RPTs establish a link between the reservoir properties (e.g., porosity, fluid saturation, clay content, etc) and the elastic properties, such as velocity, density, impedance, and wet-rock stiffness moduli (Carcione and Avseth, 2015;Hao et al., 2016;Li et al., 2019;Tan et al., 2020). Pang et al. (2019) used the P-wave quality factor and impedance to build RPTs, and shown that the estimated porosity and saturation are consistent with actual gas production results and well-log data. Gegenhuber and Pupos (2015) used RPTs based on V P /V S and acoustic impedance applied to carbonate samples.
We used V P /V S , impedance and density to build 3D RPTs, based on the EIAS model and the Zener frequency kernels, to characterize the relation between elastic properties and reservoir attributes, i.e., equant porosity and crack porosity and aspect ratio (Carcione, 2014). The templates were calibrated with ultrasonic, log and seismic data. Thereafter, a quantitative interpretation of those attributes can be performed by superposing the seismicinversion data on the templates.

LABORATORY EXPERIMENTS The Study Area
The reservoir in the Xujiahe formation of the western Sichuan Basin (China) belongs to the braided river delta sedimentary system. The main mineral components are quartz, feldspar, rock debris, mica, and heavy minerals. The size of the mineral particles is mainly medium to fine, most sandstones are sorted well, and the crushed particles are medium to poorly rounded. These sandstone reservoirs, of low porosity and permeability, have a high crack density and a heterogeneous distribution of pore fluids, due to long-term diagenesis.

Rock Specimens and Data
Fifteen samples were collected at approximately 2 km depth. They have low porosity and permeability, and are mainly composed of quartz, feldspar and rare clay. Moreover, there are calcite and quartz cementation in intergranular pores, and cracks. The ultrasonic experiments measurements (samples TSA1-TSA15) follow the experimental set-up of Guo et al. (2009), which is composed of a high pressure vessel and several units to control the temperature and the confining and pore pressures, and perform acoustic wave testing. The acquisition rate for ultrasonic waveforms (1 MHz) is 50 M/s and the time resolution is 0.02 μs.
The tests were performed at zero pore pressure and 22°C, whereas the confining pressure increases from 1 to 35 MPa (1,5,10,15,20,25,30 and 35 MPa). The differential pressure (P d ) is defined here as the difference between the confining and pore pressures. Ultrasonic compressional and shear-wave velocities (V P and V S ) at full-gas (nitrogen) were obtained by picking the first arrivals of waveforms with 1 MHz dominant frequency. These samples have low porosity and permeability, with a maximum of 13.91% and 1.37 mD, respectively. The properties are given in Table 1. Figure 1 shows a crossplot of V P and V S for varying differential pressure. As can be seen, the velocities increase with pressure, as expected. The correlation of a linear fitting between V P and V S is excellent, the R-square is 0.956, and the slope of V P vs. V S is 1.6. Figure 2 shows ultrasonic velocities as a function of differential pressure at full gas saturation. Velocities are non-linear at low differential pressure, becoming almost linear above 25 MPa, because cracks with lower aspect ratio close first and when the differential pressure exceeds approximately 25 MPa, the stiffer pores are also affected. The velocities increase faster for samples TSA4, TSA5 and TSA8, that can be due to a higher crack density, resulting in a higher sensitivity to pressure. Similar results have been observed previously by Yin et al. (2017), David et al. (2013), Pervukhina et al., 2010, Deng et al. (2015, Sun and Goldberg (1997) with different types of sandstones.

Determination of the Experimental Crack Porosity
Cracks with lower aspect ratio close first when the differential pressure increases, leading to higher velocities and decreasing crack density. We use sample TSA13 as an example to estimate crack porosity and the pressure dependency of the total, stiff (equant) and crack porosity ( Figure 3) (The former is the sum of the last two porosities). The total porosity as a function of pressure is obtained by measurement, whereas the stiff porosity by a linear extrapolation of the high-pressure trend determined between 30 and 35 MPa. The similar method was used by Yin et al. (2017) and Pervukhina et al., 2010. The crack porosity is estimated by subtracting the stiff one from the total porosity. Figure 4 shows the crack porosity as a function of the differential pressure. When the differential pressure is less than 25 MPa, the crack porosity decreases rapidly and nonlinearly. As the pressure increases, the crack porosity exhibits linear characteristics until it approaches zero at high differential pressures.

The EIAS Model
The microscopic pore structure is important because it affects the rock skeleton properties and the fluid distribution. The EIAS model can effectively address the relation between pore geometry and the elastic attributes, and it is consistent with the Hashin-Shtrikman bounds for two-phase systems (Endres and Knight 1997;Cheng et al., 2020). Assuming a background medium containing spherical pores and cracks, the resulting expressions for the effective high-frequency elastic moduli of the saturated rock are where ϕ denotes the rock total porosity, K 0 and G 0 are the bulk and shear moduli of the mineral mixture, respectively, K f is the fluid bulk modulus and β m and ζ are intermediate variables, which can be expressed as (Berryman 1995;Mavko et al., 2009;Carcione et al., 2020;Sun et al., 2020), where a is the crack aspect ratio, c is the soft-pore (crack) fraction and 1c is the stiff-pore fraction; the crack density is ρ c 3ϕc/4πa. The coefficients P 1 and Q 1 correspond to spherical pores and P 2 and Q 2 are approximations for pennyshaped cracks. The effective low-frequency moduli (when complete fluidpressure communication occurs) are where c 0 (1 − c)P 01 + cP 02 , χ 0 (1 − c)Q 01 + cQ 02 , FIGURE 1 | Crossplot of V P and V S for varying differential pressures for 15 tight-gas sandstone samples at full gas saturation. The color bar indicates differential pressure.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 646372 and P 0n P n K f 0 , (n 1, 2), Q 0n Q n K f 0 , (n 1, 2), The detailed derivation of these equations can be found in Endres and Knight (1997). The model considers the interactions between cracks and has no restrictions on the crack density. In order to generalize the model to all frequencies, we combine the EIAS model with the Zener (standard linear solid) and Kjartansson kernels, to describe the physics with two very dissimilar approaches and check the robustness of the methodology.

EIAS-Zener Model
The Zener mechanical model can be used to described the frequency dependence of dispersion and attenuation (e.g., Carcione et al., 2012;Carcione 2014). This model satisfies the Kramers-Kronig relations (Carcione et al., 2018). The minimum quality factors of the bulk and shear relaxation peaks are given by respectively, whereas the bulk and shear complex moduli are  Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 646372 where i −1 √ , f 0 is the frequency of the relaxation peak, such that for f → ∞, K → K sat HF and for f → 0, K → K sat LF . The phase velocity and quality factor of the body waves are and respectively, where v denotes v P or v S , being the frequencydependent complex P-and S-wave velocities (Carcione, 2014), respectively, where ρ is the mass density. Kjartansson (1979) presented a constant-Q model specified by the phase velocity at a reference frequency f 1 and the value of Q. To combine this model with the EIAS model, we consider f 1 50 Hz, which is the average frequency of the event corresponding to the target layer, i.e., ω 1 2πf 1 100π rad/s. The phase velocity at low frequencies is

EIAS-Kjartansson Model
For f 2 1 MHz (the ultrasonic experimental frequency), ω 2 2πf 2 2 × 10 6 π rad/s, and the phase velocity is The phase velocity as a function of frequency is (Carcione, 2014, Eq. 2.214), as an approximation for Q>>1. We obtain after some calculations, Defining, the complex modulus is (Carcione, 2014, Eq. 2.212) Therefore, the complex velocity and quality factor are respectively.

Example
The previous models are used to relate the crack properties (aspect ratio and porosity) and the frequency-dependent elastic moduli, phase velocities and quality factors. The purpose of using two dissimilar models is to verify the robustness of the approach, i.e., if the two models give similar results, we can rely on the methodology. Let us consider sample TSA13 (porosity is 13.26%). We assume that the mineral bulk and shear modulus are 39 and 36 GPa, respectively (Mavko et al., 2009), the fluid properties are obtained from Batzle and Wang (1992). We consider a Zener kernel with f 0 10 KHz and the parameters of the constant-Q kernel assumed in the previous section. Figure 5 shows the phase velocities and dissipation factors as a function of frequency, where it is clear both kernels honor the velocity dispersion between the seismic and laboratory frequency bands. Therefore, they can be used to build multi-scale 3D RPTs in the ultrasonic, logging and seismic frequency bands. Thereafter, we assume a crack porosity of 0.2% and crack aspect ratios of 0.001, 0.0014, 0.0018, 0.0022, 0.0026. Velocities and dissipation factors as a function of frequency are shown in Figures 6A,B. When the crack porosity is a constant and the crack aspect ratio increases, the rock becomes stiffer and the P-wave velocity increases, but attenuation decreases. If we take a crack aspect ratio of 0.002 and crack porosities of 0.2, 0.3, 0.4, 0.5, 0.6%, the results are shown in Figures 6C,D, where the trend opposite to the previous one.

Building the RPTs
By setting the crack aspect ratio, crack porosity and total porosity as variables, the corresponding P-wave impedance, density and V P /V S value are calculated by using the EIAS-Zener model. Then, we build the 3D RPTs at ultrasonic (1 MHz) well-log (10 kHz) and seismic (50 Hz) frequencies, and the seismic template is used to estimate the microstructure of the sandstone reservoirs. Similar templates have been used by Pang et al. (2020a) to estimate the microfracture porosity of deep carbonate reservoirs. Figures  7A,B show the 3D RPTs at ultrasonic and seismic frequencies, respectively (the well-log template is similar to the seismic one).
As the stiff porosity increases, the P-wave impedance and mass density decrease, but the variation in V P /V S is small. With the increase of crack porosity, the P-wave impedance decreases and V P /V S increases, but the density hardly changes. As the crack aspect ratio increases, the P-wave impedance increases, V P /V S decreases. This behavior is similar at all frequencies.

Calibration at Ultrasonic Frequencies
The templates are calibrated with the crack porosity obtained in Section Determination of the Experimental Crack Porosity. Figure 8 superposes the data to the template. The porosity (2.43-13.91%) of samples is basically consistent with the corresponding values in the templates, and the agreement is also good for the crack porosity. What's more, the average error between the predicted porosities (based on ultrasonic wave attributes) and the measured one is 0.00455. In addition, as this increases, the P-wave impedance decreases. Figure 9 shows the calibration of the 3D RPTs at log and seismic frequencies, where the data is extracted from well B. The template at seismic frequencies is calibrated with inversion results (P-wave impedance, density and V P /V S ) of pre-stack seismic data that has been verified, whole the porosity is taken from the logs. Similar methods have been used by Pang et al. (2020b). The velocities and density at seismic frequency band are extracted from the seismic traces around the borehole by using a three-term inversion method (Aki and Richards, 1980). We use the partial-stack method to obtain the seismic data based on pre-stack angle gathers, which improves the signalnoize ratio and suppresses random noises. In addition, the velocities and porosity at sonic-log frequency band are extracted from the well-log data by using the velocity-time semblance method, which avoids picking arrival times and can accurately determine the phase velocities for different modes (Wang et al., 2020), and the density is obtained from density logging. Figure 9 shows that the porosity is basically consistent with the corresponding values in the templates. Figure 10 shows a crossplot of V P and V S of well B, where the velocities increase with decreasing porosity, as expected. The Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 646372 6 correlation of a linear fitting between V P and V S is excellent, and the slope of V P vs. V S is 1.12.

Calibration at Well-Log and Seismic Frequencies
The P-wave impedance, density and V P /V S decreases with increasing porosity, and the latter decreases with increasing crack aspect ratio. The porosity is low (less than 7%), which is consistent with the geological characteristics of the reservoirs.

ESTIMATION OF THE MICROSTRUCTURE FROM SEISMIC DATA
The predictions of reservoir stiff porosity and crack porosity and density are performed on a 2D seismic line crossing wells A and B. P-wave impedance, density and V P /V S are overlapped onto the 3D RPTs, and the microstructural properties are estimated. Figure 11A shows the porosity from well log A and B are superposed, showing a good agreement. By taking well B as an example, the average error between the predicted porosities (based on seismic wave attributes) and the measured porosities is 0.00923, which is higher than the one at ultrasonic frequency, because the resolution of the seismic data is low. Figures 11B,C show the crack porosity and density, respectively. It shows that well B has a higher crack porosity and density. The horizontal sections are created from target layer (depth 4,787 m) of the study area in the Xujiahe formation of the western Sichuan Basin (China). Figures 12, 13 show horizontal section of those attributes, which are simulated by using EIAS-Zener and EIAS-Kjartansson models, respectively. The results show that the two models give similar results, thus indicating that the methodology is working. The total porosity and crack porosity and density effectively serve to identify highquality reservoirs. Well B is located in an area with high porosity and crack content, where the target layer has a good pore connectivity, storage capability, and natural gas accumulation conditions. The gas production report of well B indicates 5148929 m 3 per day, while well A yields 57559 m 3 per day. The results agree with the actual production of the two wells. Therefore, the methodology is reliable, and is expected to provide an important basis for further rock-physics studies and field applications of seismic exploration.

CONCLUSION
Rock-physics templates constitute an easy-to-use "tool-box" for pore microstructure properties interpretation from seismic data. We have implemented the EIAS-Zener model to describe the crack properties and elastic moduli of tight-gas sandstone reservoirs as a function of frequency from the seismic to the ultrasonic bands. Since P-wave impedance, mass density and V P /V S are sensitive to the pore and crack properties, we built the templates to estimate these attributes from seismic data. First, we calibrate the templates with ultrasonic, well-log and seismic data. Then, we estimate the total porosity and crack porosity and density from inversion of seismic data. The results indicate that these estimated attributes are consistent with the actual gas production reports. Thus, these templates can effectively be used to estimate the microstructure characteristics of tight-gas sandstone reservoirs based on seismic data.

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

AUTHOR CONTRIBUTIONS
WC: Methodology, writing original manuscript JB: Funding support, data and discussion, writing original manuscript, supervision JC: Model building, revising manuscript MP: Seismic data analysis CW: Ultrasonic experimental data analysis.