An Efficient Rock Physics Scheme for Estimating Crack Density and Fluid Saturation of Shale Gas Reservoir

We propose a simple rock physics model for the characterization of elastic properties of shale. The model combines a dual-porosity concept and the effective medium theory for constructing the anisotropic elastic tensor of the multimineral organic-rich shale. Based on the model, we address how to estimate two key shale gas evaluation parameters, i.e., crack density and gas saturation from well-log and seismic data. Application to Wufeng-Longmaxi Shale shows that the estimated crack porosity decreases with increasing burial depth and decreasing clay content. The analysis indicates that the microcracks are mainly developed among clay minerals, which is consistent with the results from mercury injection and SEM imaging experiments. More importantly, we show that the velocity of the Wufeng-Longmaxi Shale is primarily controlled by the crack porosity instead of the total porosity. Both P- and S-wave velocities decrease linearly as the volume of microcrack increases. The fluid substitution analysis shows that the Poisson’s ratio and P-impedance of the shale are sensitive to the change of pore-fluid saturation. Based on the above sensitivity analyses, we customize a rock physics template for quantifying crack density, and gas saturation from the shale elastic properties. The interpretation results show that there is an overall good agreement between the measured and predicted petrophysical properties from well-log and seismic data.


INTRODUCTION
Quantitative seismic interpretation is the core of shale gas reservoir identification and evaluation. During the exploration stage, seismic data are employed to determine buried depth, thickness, and structural shape of the shale reservoir (Rich and Ammerman, 2010;Bachrach et al., 2014;Jiang et al., 2019). This is followed by applying the inversion technique for mapping the distribution of organic matter content, porosity and this leads to the determination of the favorable area with economical gas potential (Chopra et al., 2012). In the development stage, the seismic technology is used to determine the mechanical and anisotropy characteristics of the reservoir, so as to provide the basis for the deployment of horizontal wells, and well bore design and fracturing stimulation (Goodway et al., 2010;Sena et al., 2011). Understanding the relationship between the elastic and petrophysical properties of shale is fundamental to the geophysical research of shale gas exploration and development. Shale gas reservoir has complex lithologic, pore and structural characteristics. Among them, the compactness caused by low porosity and low permeability is a more obvious feature of shale gas reservoir (Passey et al., 2010). Former empirical formulas and theoretical rock physics models developed for conventional clastic and carbonate formations may not be applicable for shale gas reservoirs.
A lot of work have been done around the experimental and theoretical investigations of shale elastic properties. Vernik and Nur (1992); Vernik and Liu (1997) measured the elastic anisotropy of Bakken Shale under dry conditions. They analyzed the effects of organic matter content and maturity on rock velocity and anisotropy. They pointed out that the kerogen has a dramatic effect on seismic velocity and the anisotropy of shale is largely controlled by the preferred direction of minerals and micro-fractures. Sondergeld et al. (2000); Sondergeld and Rai (2011) conducted experimental study on the acoustic properties of Kimmeridge Shale and found that the anisotropy of shale increases with increasing organic matter content. The increase of organic matter content will lead to the decrease of density, resulting in the effect opposite to compaction. They also showed that the assumption of weak anisotropy (Thomsen, 1986) cannot be used in the seismic modeling of shale. Domnesteanu et al. (2002); Deng et al. (2009) studied the velocity and attenuation of shale under different saturation and pressure conditions using ultrasonic transmission technique. The experimental results suggest that there is a strong link between the pressure-dependent anisotropy of velocity/attenuation and the pore geometry, connectivity and the response of pore fluid to the propagation direction of the seismic waves. Hornby et al. (1994) combined the self-consistent approximation (SCA) and the differential effective medium (DEM) theories to characterize the seismic elastic properties of Kimmeridge Shale. In the model, clay is treated as the supporting matrix that forms the continuous skeleton of the rock, other minerals such as quartz, feldspar, and pyrite are regarded as isolated inclusions dispersedly distributed in clay. Carcione et al. (2011) used Backus average to evaluate the effect of clay particle orientation on the elastic properties and anisotropy of shale. Guo et al. (2013) developed a seismic petrophysical model of organicrich shale and investigated the effects of pore shape and mineral composition on the elastic properties of shale.
Organic-rich shale exhibits relatively higher resistivity, higher sonic transit-time and lower density than nonorganic rich sediments (Passey et al., 2010). Tight shale hydrocarbon reservoirs are characterized by total porosity in the range of 3-7% and matrix permeability in the range of micro-to nano-Darcy. The seismic elastic properties of shale are controlled by its own structural characteristics. The key to analyze its seismic elastic properties and influencing factors is to accurately describe the structural characteristics of shale at microscale including the spatial distribution of main constituent minerals, organic matter, and pore structure (Vasin et al., 2013;Luan et al., 2016). Due to the differences in sedimentary history and depositional environment (i.e., stress field change, mineral composition, etc.) of different reservoir rocks, the petrophysical experimental results for specific reservoirs are also regional, and can not be extrapolated arbitrarily Deng et al., 2021). Therefore, it is important to establish site-specific rock physics model to characterize the local elastic properties of shale.
Natural cracks in shale provide necessary fluid conduits for shale gas migration and can greatly enhance the seepage capacity (Guo et al., 2014;Wu et al., 2020). The degree of microcracks contributes to the total porosity which affects the storage capacity of shale reservoir. The evaluation of microcrack volume is of great importance for the exploration and development of shale gas. On the other hand, shale gas saturation is important for predicting reservoir performance in terms of production (Lucier et al., 2011;Qi et al., 2017). Accurate estimation of gas saturation remains as one of the most difficult tasks associated with shale reservoir characterization (Rezaee, 2015;Qi et al., 2021). In order to predict the rock and fluid parameters of gas shale by geophysical methods, it is necessary to establish the relationship between elastic properties and petrophysical parameters. In this work, we address an important issue of estimating two key shale gas evaluation parameters, i.e., crack density and gas saturation from well-log and seismic data. The paper is organized as follows: first, we introduce a simple but effective rock physics modeling scheme for organic-rich shale. Then, we demonstrate how this model can be applied to the inversion of pore structures, fluid substitution, and reservoir characterization with rock physics template. This is followed by the discussion on the novelty and applicability of the model. Finally, we elaborate the main conclusions and set the basis for future works.

METHODOLOGY
In this work, we propose a simple anisotropic dual-porosity effective medium model (ADPM) for calculating the elastic property of the organic rich shales. The microstructure of the organic-rich shale can be extremely complex (Vanorio et al., 2008;Uvarova et al., 2014;Qi et al., 2021). The pores of different shapes, such as intergranular pores and microcracks, can develop between or on the edges of the rock matrix minerals such as quartz, and clay. Intraparticle pores of nanoscale can also develop among the organic matter. Despite the complex nature of the shale microstructure, we assume that at any specific logged depth location, the shale porosity system (ϕ total ) can be simplified and divided into two spaces: the stiff-pore space (stiff-porosity) and the microcrack space (crack-porosity), i.e., ϕ total ϕ stiff + ϕ crack . (1) Stiff porosity (ϕ stiff ) and crack porosity (ϕ crack ) are defined here as the volumetric fraction of spheres with larger aspect ratio and oblate spheroid inclusions with very small aspect ratio, respectively. The key to our assumption is that we do not consider how the different pore structures fit into the rock frame since their distributions can be arbitrary and are unknown. Later we show that this assumption greatly simplifies the modeling (inversion) procedure and yields results that are consistent with laboratory and field measurements. With the above dual-porosity model, we propose the following rock physics workflow for modeling shale elastic properties. First, since different components in the organic-rich shale possess very different stiffness, an accurate mineralogy analysis is crucial for constructing the right rock matrix to start with. First of all, we can obtain the volume fraction of the mineral content in shale from the spectral gamma-ray (well-log) or the X-ray diffraction analysis (core sample). When the exact mineralogy of shale cannot be assessed, simple gamma-ray log can be used to provide an estimate for the volumes of silt and clay in the formation. The mineralogy of the organic-rich shale can be broadly divided into the brittle minerals, clay and kerogen components. The brittle minerals include quartz, carbonate and pyrite whereas the kerogen is the typical soft material. When the volume fractions of each component are determined, we can apply the Voigt-Reuss-Hill average (Mavko et al., 2009) to compute the isotropic effective moduli of the brittle mineral aggregate. Because the shale anisotropy is mainly caused by the lamination of clay and kerogen, we use a three-phase Backus average (Backus, 1962;Vernik and Nur, 1992) to mix the brittle minerals, clay and kerogen. This procedure yields the anisotropic effective moduli of multimineral solid matrix.
The organic-rich shale is typically of low porosity and low permeability with poor pore connectivity. These shales may contain randomly distributed microcracks with crack size below well-log resolution. Theoretical estimates of the effective moduli of the cracked porous shale not only depend on the volume fraction of the inclusions (i.e., cracks and pores) but also the geometric details of shapes and spatial distributions of these inclusions. Several effective medium theories can be used to estimate the effective moduli of the shale. In this work, we use the anisotropic self-consistent approximation (SCA) theory to compute the stiffness tensor of the shale frame. The choice of the SCA model is due to the following consideration. Since the inclusion adding process in the SCA theory is a symmetrical process, it does not identify a specific phase as the host material. This is different for the differential effective medium (DEM) theory. Due to the difference in depositional and diagenetic processes, shale can transits from clay-supporting matrix to quartz-supporting matrix as the burial depth changes Deng et al., 2021). The reservoir section usually contains shales of more than one type. Therefore, we consider the SCA model is more flexible and can adapt to more complex lithologic enviroment. By incorporating the dual-porosity concept into the SCA modeling scheme, the effective stiffness tensor of the shale can be written as (2) where, and i represents the matrix (m) or the pore spaces (stiff, crack), I is the identity tensor. C i is the stiffness tensor matrix of the i-th phase. G i is a fourth-rank tensor calculated from the response of an unbounded matrix of the effective medium Eshelby (1957); Mura (1982). Brown and Korringa (1975) extended the Gassmann's equations and established the model suitable for fluid replacement in anisotropic media such as shale. On the basis of dry shale model, we apply the Brown-Korringa's (BK) equation to compute the stiffness tensor of the saturated shale. In the BK model, the compliance of the saturated shale is given by where S dry ijkl is the effective elastic compliance tensor element of dry rock, S sat ijkl is the effective elastic compliance element of rock saturated with pore fluid, S 0 ijkl is the effective elastic compliance element of the solid mineral, β fl is the fluid compressibility, and β 0 is the mineral compressibility. The stiffness can be obtained by taking the inverse of the compliance matrix. To model the partially saturated gas shale, we use the Brie's model (Brie et al., 1995) to calculate the effective fluid moduli of a brine and gas mixture. In accordance with the BK model, the compliance of the fluid mixture can be calculated as The fluid compliance depends on the bulk moduli of brine K w and gas K g , the gas saturation S g and an empirical exponent e. The empirical exponent e varies from 1 for patchy saturation to 3.4 for uniform saturation. They provide the upper and lower bounds for the effective bulk modulus of the mixed fluids, respectively. We can substitute the equivalent compliance of brine-gas mixture (5) into the BK model (4) to replace the empty pores with fluids.
Based on the above anisotropic dual-porosity model, the equivalent elastic tensor of fluid saturated rock can be converted to velocities according to the relationship provided by Thomsen (1986). For a vertical well, the acoustic logging tool records the P-and S-waves traveling perpendicular to the bedding plane, that is parallel to the axis of symmetry of the transversely isotropic shale formations. The expressions for these velocities are where C 33 and C 44 are the stiffnesses in Voigt notation. ρ is the effective density of the shale, which can be calculated by the volume average of the grain and fluid densities To test the effectiveness and applicability of the APDM model, in the following sections, we apply the rock physics model to a suite of analyses including pore type inversion, shear-log prediction, fluid substitution, and seismic reservoir characterization.

Geological Background
We conduct rock physics analysis on the organic-rich shale in the Wufeng-Longmaxi Formation. The Wufeng-Longmaxi Shale are widespread units, and their geologic characteristics are similar to five of the leading shale gas producing basins in North America (Pollastro, 2007;Guo and Liu, 2018;Xu et al., 2019). Therefore, establishing a quantitative characterization model for the Wufeng-Longmaxi Shale would have global impacts. Well-Y1 is a shale gas exploration well located in a shale gas field in Southern Sichuan Basin, China. The section in Figure 1 is composed of the lower Silurian Longmaxi formation and the upper Ordovician Wufeng formation. The organic-rich shale at the lower part, which shows high gamma and low density, and is the main target layer. The thickness of organic-rich shale in well-Y1 is 67 m, the total organic carbon content (TOC) is 1.8-4.8 wt %, and the thermal maturity is 2.5-3.0. All geological evaluation parameters show that Wufeng-Longmaxi Formation in the work area has superior conditions for shale gas accumulations. As shown in Figure 2, measurements of density, gamma ray, shallow-and deep-resistivities, compressional, and shear velocities are available. The porosity can be calculated based on the density log. In these reservoirs, the total porosity varies from 2 to 10% with an average value of 5%. Figure 2A shows that the average values of clay, quartz, and carbonate are 38, 48, and 11%, respectively. Pyrite is also present, ranging from trace amounts to as much as 7%. The quartz is present in the form of silt and siliceous cement. TOC is estimated from resistivity and sonic logs using the method of Passey et al. (1990). The Wufeng-Longmaxi shale reservoir is a typical low-resistivity reservoir as indicated by Figure 1E. Conventional resistivity method, such as Archie's equation (Han T.-C. et al., 2021), will lead to overestimation on the water saturation as shown in Figure 1F.
To overcome this problem, water saturation is estimated using the neutron-density porosity difference method (Qi et al., 2021). Figure 2D shows that the corrected water saturation value varies from 18 to 100% in the reservior section. In the rock physics modeling, the bulk moduli and densities of the brine and gas of in-situ pressure and temperature are estimated using the equations of Batzle and Wang (1992). The elastic moduli and densities for the minerals and pore-fluids are presented in Table 1.

Inversion of Different Pore Volume Fractions
The rock physics model provides a nonlinear relation between Pand S-wave velocities and the volume fractions of soft and stiff pores. Using this relation, we can forward calculate the P-and S-wave velocities of rocks based on rock and fluid properties including mineral moduli, porosity, and fluid saturation and pore aspect ratio. Alternatively, we can determine the stiff and soft porosities by matching well-log velocity data with theoretically predicted velocities. Here, we assume that the stiff pores and microcracks are spheroidal inclusions with aspect ratios of 1 and 0.01, respectively. Thus, we keep the pore aspect ratio as constant at all depths in the inversion procedure. The choice of aspect ratio 1 for stiff pore is based on the fact that sphere is the stiffest among all pore shapes. For cracks of small aspect ratio, the key parameter is the crack density ϵ 3ϕ crack /4πα crack (O'Connell and Budiansky, 1974), which is proportional to the crack porosity divided by the aspect ratio. The choice of exact crack aspect ratio value is of secondary importance. We also assume that all cracks have identical shapes, i.e., aspect ratios. In this case, the effect of crack porosity on elastic moduli is equivalent to the crack density (Vernik, 2016). To determine the volume fractions of cracks and stiff pores from sonic-log data, we propose the following objective function where V obs and V model are the measured and modelled velocities, respectively. W p , W s are the weights, for case that the shear-log is not available, W p 1, W s 0 and when both compressionaland shear-logs are available, W p 0.5, W s 0.5. By minimizing the objective function (9), ϕ crack and ϕ stiff are jointly inverted from two velocities V p and V s . A grid search method can be applied to find the two porosities required to accomplish the velocity match at each well-logging depth.
Other formation properties such as mineral component, total porosity and fluid saturation are provided by the welllog data.
The estimated crack and stiff porosities are plotted in Figure 2G. We see that the microcracks are well developed in the upper and middle parts of the Wufeng-Longmaxi Shale. As the burial depth increases, the amount of microcracks diminishes and the pore-space is dominated by stiff pores. The average porosities of the crack and stiff pores in the reservoir section are 0.012 and 0.043, respectively. The volume fraction of the stiff pores is 3.6 times of the microcrack overall. Figures 2E,F shows the modeled P-and S-wave velocities based on the inverted crack and stiff porosities. There is a satisfactory match between the measured and the modeled velocity. The root-mean-square error (RMSE) for the predicted P-and S-wave velocities are 34.9 m/s and 23.4 m/s, respectively. The correlation coefficients between the measured and predicted velocities are 0.992 3 for P-wave and 0.991 6 for S-wave. The  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 829244 consistent prediction of the two velocity logs shows that the rock physics model is effective for characterizing the elastic behavior of the shale. To further analyze the relationship between the elastic properties of shale and its pore structure, we plot the velocities as functions of different porosities. Figure 3A shows that there is a weak dependence of the velocities on the total porosity of the shale. The P-wave velocity increases slowly with increasing total porosity, which is counter-intuitive. The total porosity shows limited influence on the velocities. Figure 3B shows the relationship between the velocities and the crack porosity. Interestingly, both P-and S-wave velocities exhibit overall linear correlations with the crack porosity. A linear leastsquare fit provides the following relations The velocities decrease as the volume of microcracks in shale increases. The clear dependence of the velocity on crack porosity shows that the presence of small-aspect-ratio cracks is a key parameter in controlling the elastic property of the Wufeng-Longmaxi Shale. Figure 4 shows that the development of the microcracks in the Wufeng-Longmaxi Shale generally decreases as burial depth increases. The data points are color-coded by the clay content, as can be seen the crack porosity tends to increase with increasing volume of clay in shale. This observation is consistent with the laboratory results of Sun et al. (2019) who conducted mercury injection experiments to quantify the volume of microcracks in the Longmaxi shale samples. They found that the majority of the microcracks were formed inside the clay particles, followed by those developed in organic matter and brittle minerals (see more details in the discussion section). The SEM images in Figure 5 reveal the microstructure of the Wufeng-Longmaxi organic-rich shale. In Figures 5A,B, the pores between clay platelets or around the edge of rigid particles are shown as elongated, compliant shape. Whereas the intergranular pores, which have a rigid triangular structure, are mainly developed between the quartz minerals as shown in Figures 5C,D. Our observation that the microcracks are more likely to develop in the clay is consistent with both the mercury injection statistics and the actual microstructure revealed by the SEM images of the Wufeng-Longmaxi Shale. Once we determine the crack and stiff porosities of the shale formation, a detailed fluid substitution analysis can be carried out at each logged depth. This is discussed in the next section.

Fluid Substitution Analysis
Fluid substitution is important for examining the sensitivity of the elastic properties to the change of amount of fluids in the porespace. It directly helps seismic interpreter to determine area with favorable gas potential from cross-plot of seismic attributes. With the calibrated crack and stiff porosity logs, we can conduct fluid substitution for the entire reservoir section based on the dualporosity model. By assigning different levels of gas saturation S g to Eq. 5 and substituting the effective fluid compressibility into Brown-Korringa's equation, we can evaluate the elastic properties of the partially saturated shale. It has been confirmed by studies that due to low permeability and poor pore connectivity, patchy saturation instead of uniform saturation is likely to take place in tight reservoir rocks (Si et al., 2016;Li et al., 2020). We pick the exponent e 1 in Brie's model based on an empirical fit to the data. This corresponds to the upper bound for the velocitysaturation relation in shale.  The water saturation log in Figure 2D indicates that the gas saturation in the shale formation varies from 0 to 72% with an average value of 37%. As shown in Figure 6A, the corresponding Poisson's ratio varies from 0.25 to 0.36 and the P-impedance varies from 0.95e 4 -1.2e 4 g ·m/cm 3 ·s. The red squares in Figure 6 represent the shale formation with full water saturation. The blue circles in Figures 6A-D denote the original well-log data and the formation with gas saturation of 30, 70 and 90%, respectively. For low formation gas saturation, the gas shale shows similar characteristics with the water shale and it is difficult to distinguish them on a Poisson's ratio versus P-impendance cross-plot. When the gas saturation increases to 70%, both the Poisson's ratio, and P-impedance decrease dramatically. For 90% gas saturation, the Poisson's ratio becomes less than 0.25 and the P-impedance reduces by 9%. Overall, with favorable gas saturation larger than 70%, there exists a clear separation between the original, and substituted results on the cross-plot. This indicates that the Wufeng-Longmaxi shale is relatively sensitive to the change of pore fluid saturation. The results suggest that the cross-plot between Poisson's ratio and P-impedance can be useful and applied to seismic inversion results for gas zone identification in the work area.

Seismic Characterization With Rock Physics Template
Rock physics template (RPT) is an essential tool for mapping the effect of fluids, porosity, and lithology on elastic properties such as P-, S-impedances and V p /V s (Avseth et al., 2010;Zhou et al., 2021). Because impedance and V p /V s are properties that can be obtained from seismic inversion, RPT allows quantitative interpretation and classification not only of well-logs, but also of seismic data. The main goal is to use the RPT for quantifying two primary properties, i.e., crack density and gas saturation from the well-log and seismic data in the work area. In the last two sections, we have shown that the velocities of the Wufeng-Longmaxi Shale strongly linked to the crack porosity while the Poisson's ratio is sensitive to the fluid saturation. The P-impedance and V p /V s can be viewed as proxies for P-wave velocity and Poisson's ratio, respectively. Therefore, we design the RPT so that it outputs a cross-plot between the P-impedance and V p /Vs for quantitative interpretation. To construct the RPT, we set the crack density and gas saturation as variables in the model and fix other depth-dependent properties to reference values. For example, model parameters such as effective mineral moduli, porosity and fluid properties are defined using the average values in the reservoir section. Thus, for each crack density and gas  : (A,B) shows that the pores between clay platelets or around the edge of rigid particles are of elongated, compliant shape; (C,D) shows that the intergranular pores, which have a rigid triangular structure, are mainly developed between the quartz minerals.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 829244 saturation, the template predicts a pair of P-impedance and V p /V s . Figure 7 shows the RPT with the data points color-coded by the well-log crack density. The RPT predicts that at high gas saturation, the V p /V s decreases with increasing crack density.
When the gas saturation is below 60%, the V p /V s increases as crack density increases. The P-impedance decreases with increasing crack density at all saturations. There is a good agreement between the data and the predictions by the RPT. FIGURE 7 | Rock physics template analysis for the interpretation of crack density and gas saturation from cross-plot between the V p /V s and P-impedance. The data points are from Well-1 and color-coded by crack density log. The arrow shows the direction of increasing crack density.
FIGURE 8 | Rock physics template analysis for the interpretation of crack density and gas saturation from cross-plot between the V p /V s and P-impedance. The data points are from Well-1 and color-coded by gas saturation log. The arrow shows the direction of increasing gas saturation.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 829244 8 FIGURE 9 | Seismic inversion results of the shale gas field: (A) acoustic impedance and (B) shear impedance. The location of Well-Y1 is indicated by the black line.
FIGURE 10 | By projecting the inversion results in Figure 9 onto the rock physics template in Figures 7, 8, we obtain the interpreted sections of (A) crack density, and (B) gas saturation.
Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 829244 9 The crack density of the shale formation can be accurately predicted from the RPT analysis. Figure 8 presents the same RPT with the data points color-coded by the well-log saturation. Again, there is an overall excellent match between the log and predicted saturation values. The interpretation of the well-log data shows the effectiveness and accuracy of the constructed RPT. Next, we apply the RPT to the seismic data and obtain the spatial distribution of the crack density and gas saturation.
To obtain the seismic attributes required for RPT analysis, we perform the pre-stack seismic inversion. In the inversion process, the differences between the seismic amplitudes from different reflection angles are used to obtain the corresponding variation in the shale elastic properties. The pre-stack inversion is based on the Fatti's reflectivity equation in terms of acoustic and shear impedances (Fatti et al., 1994). The inversion results are shown in Figure 9. The top and bottom of the shale reservoir are denoted by the horizons T1 and T2, respectively. The highquality shale at the bottom of Longmaxi formation in the work area are characterized by low impedance (blue color stripes). We can compute the V p /V s by taking the ratio of Pand S-impedance volumes. Similar to the log data analysis, by projecting the seismically-derived V p /V s , and P-impedance onto the RPT, we can predict the corresponding crack density and gas saturation. Figure 10 shows the interpreted sections. The zones of interest can be clearly identified from the interpreted sections. The main target at the lower part of the Wufeng-Longmaxi Formation shows high gas saturation and considerable microcrack concentration. The maps indicate that the lateral distribution of high-quality shale is stable and distributed throughout the region. Such depositional pattern is consistent with the geological studies of this area (Chao et al., 2012;Liu et al., 2013). The difference in the lateral change between crack density and gas saturation can be explained by the fact that the rock and fluid property variations contribute to the change in P-and S-impedances differently (Avseth et al., 2010;Vernik, 2016). The study shows that the RPT based on dualporosity model is an efficient tool for quantitative interpretation of the well-log and seismic data of the Wufeng-Longmaxi Shale reservoir.

DISCUSSION
Based on the dual-porosity model, we proposed an inversion scheme for estimating the volume fractions of different pore structures. Our analysis shows that the majority of the microcracks are formed within the clay minerals, which is consistent with the observation from the SEM images (see Figure 5). Sun et al. (2019) carried out a thorough analysis on the correlation between the volume of microcrack and the mineral components in the Longmaxi shale samples based on the electron microscopy and high-pressure mercury injection experiments. The results show that the contribution of clay minerals to the volume of microcracks is the largest, accounting for 70%, followed by the organic matter, and accounting for 24%. The brittle minerals (including quartz, feldspar, and dolomite) contribute to the least volumes of microcracks, accounting for 6%. The experiment also shows that there exists a positive correlation between the volumes of microcrack and clay. This agrees with the sonic-log inversion results based on the dual-porosity model.
The clay minerals reduce the brittleness of shale to a certain extent, this is not in favor of forming the microcracks by mechanical deformation due to tectonic stress (Han et al., 2021b). However, the microcracks can form in terms of small openings between the crystal layers of different clay minerals as revealed from the SEM images . On the other hand, clays are composed of fine sheet-like particles, they tend to form pores with much smaller aspect ratios than those associated with quartz grains. These pores are very compliant, most of which can be closed by the lithostatic stress at larger burial depth. However, overpressure is a common phenomenon in the shale gas reservoirs of Upper Ordovician Wufeng Formation and Lower Silurian Longmaxi Formation in the southern Sichuan Basin (Gao et al., 2019;Han et al., 2021a). Hydrocarbon generates from organic matter during thermal evolution, which is the main cause for abnormal high pressure in the region. The overpressure reduces the effective stress on reservoir and therefore, the majority of the microcracks still remain open for the Wufeng-Longmaxi Shale. This corroborates with our observation of considerable microcrack growth in well-Y1. Xu and White (1995) published a velocity model for clay-sand mixtures, which is developed in terms of the Kuster-Toksöz effective medium theory and Gassmann model. In their model, the total pore-space is assumed to consist of pores associated with sand grains, and clays, respectively. The model also assumes that the pores of small aspect ratio are strictly developed in the clays whereas the pores of large aspect ratio only exist among the sand grains. Therefore, for its application one must first determine the volume of sand grain and clay content. In presence of more complex mineral components and microstructures, such as in organic-rich shales, and the classic Wu-White's model may no longer be applicable. Unlike the Xu-White's model, the dualporosity model proposed in this work does not differentiate how the pore structures fit into the rock frame. Instead, we consider that the microcrack and stiff pores can develop in all mineral components. This assumption greatly simplifies the modeling and inversion procedures. The results as discussed earlier confirms the effectiveness of the dual-porosity model. The pore-volume inversion workflow can serve as fast evaluation tool for assessing the microcrack development in shale reservoir. In this work, the model is applied to the Wufeng-Longmaxi Shale Formation which is of uniform lithology. More future work is needed to examine the applicability of the model for reservoirs with complex lithologies.

CONCLUSION
In this work, we propose a simple dual-porosity model for the characterization of elastic properties of shale. The model assumes that the pore-space of shale can be divided into two parts including the pore-spaces of microcracks and stiff pores, respectively. The crack and stiff pore are defined as spheroid Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 9 | Article 829244 inclusion with contrasting aspect ratios. The dual-porosity concept is incorporated into the anisotropic self-consistent approximation theory for obtaining the effective elastic tensor of the multimineral organic-rich shale. We use the model to analyze the elastic behavior of the Wufeng-Longmaxi Shale in Sichuan Basin, China. We address the important issue of estimating two key shale gas evaluation parameters, i.e., crack density and gas saturation from well-log and seismic data. The main conclusions can be summarized as follows: 1) The log data analysis shows that the estimated crack porosity in Wufeng-Longmaxi Shale exhibits a decreasing trend with increasing burial depths and decreasing clay content. The microcracks are mainly developed between clay minerals, which is consistent with the experimental results and observation from SEM images.
2) The velocities of the Wufeng-Longmaxi Formation is primarily controlled by the crack porosity instead of the total porosity. Both the P-and S-wave velocities show linear decreasing trends with increasing crack porosity in shale.
3) The fluid substitution analysis indicates that the gas shale is sensitive to the change of pore-fluid saturation. We can separate the gas shale from the water shale on a cross-plot between the Poisson's ratio and P-impedance.
4) The rock physics template predicts that at low gas saturation, V p /V s increases with increasing crack density whereas at high gas saturation, V p /V s decreases with increasing crack density. The well-log crack density and gas saturation can be accurately predicted from the RPT analysis. The interpreted seismic sections show that the high-quality shale at the bottom of Wufeng-Longmaxi Formation is stable and distributed throughout the region.

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.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for submission.