The Lower Silurian shale-gas formation in the south of the Sichuan Basin represents a strong transverse isotropy with vertical axis of symmetry (VTI) feature. Successful characterization of shale-gas formation requires handling the great influence of anisotropy in the seismic wave propagation. Seismic amplitude variation with offset (AVO) inversion for VTI media using PP-waves only is a difficult issue because more than three parameters need to be estimated and such an inverse problem is highly ill posed. We have applied an AVO inversion method for VTI media based on a modified approximation of the PP-wave reflection coefficient. This approximation consists of only three model parameters: the acoustic impedance (attribute A), shear modulus proportional to the anellipticity parameter (attribute B), and the approximated horizontal P-wave velocity (attribute C), which can be well-inverted and have great interpretation capability in shale-gas reservoir characterization. A statistical-rock-physics method was then applied to the inverted attributes for quantitative interpretation of the shale-gas reservoir. A Markov random field is combined with Bayesian rule to improve the continuity and accuracy of the interpretation results. Shales can be successfully discriminated from surrounding formations by using the attribute pair A-C, and the organic-rich gas-bearing shale can be successfully identified by using the attribute pair C-B. Comparison between the prediction results and well logs demonstrates the feasibility of the inversion and quantitative interpretation approaches.


Shales represent intrinsic elastic anisotropy (transverse isotropy with vertical axis of symmetry [VTI]) that is much stronger than other sedimentary rocks. Vernik and Nur (1992) find that shale anisotropy is attributed to the presence of organic matter. Hornby et al. (1994), Lonardelli et al. (2007), and Wenk et al. (2007) suggest that such intrinsic anisotropy is mainly caused by preferentially orientated clay platelets and cracks. Anisotropy has a significant influence in seismic wave propagation. Wright (1987), Kim et al. (1993), and Thomsen (1993) analyze the effects of anisotropy on seismic amplitudes for VTI media. Keith and Crampin (1977), Daley and Hron (1977), Graebner (1992), and Schoenberg and Protázio (1992) determine analytical solutions of reflection coefficients for VTI media, and varied forms of approximations have been proposed in many investigations, e.g., Banik (1987), Thomsen (1993), Ursin and Haugen (1996), Rüger (1997), Zillmer et al. (1998), Vavryčuk and Pšenčík (1998), Stovas and Ursin (2003), Shaw and Sen (2004), and Zhang and Li (2013). However, the amplitude variation with offset (AVO) inversion for VTI media using PP-wave only is a difficult issue. Plessix and Jonathan (2000) study the feasibility of anisotropic parameter estimation from reflection coefficient of VTI media and demonstrate that the recovery of all five parameters is difficult in practice. Stovas et al. (2006) invert AVO intercept and gradient for net-to-gross and oil-saturation estimation for VTI media. Lin and Thomsen (2013) present a method to extract the anisotropy parameter δ based on the difference between real (measured) and synthetic seismic amplitudes.

The purpose of seismic reservoir characterization is essentially to identify the sweet spots that represent the most valuable target drilling areas. True seismic amplitudes are diagnostic of lithology and hydrocarbon anomalies. Mukerji et al. (2001) and González (2006) propose quantitative interpretation approaches for seismic reservoir characterization by using seismic amplitude and related attributes. Identifying sweet spots of shale oil and gas includes brittleness and organic-matter content estimation. Chopra et al. (2018) and Convers et al. (2017) estimate brittleness from inverted elastic parameters by empirical well-log analysis. Wang et al. (2019) propose a novel statistical rock-physics method for quantitative interpretation of brittleness. Chopra et al. (2018) and Kumar et al. (2018) estimate organic-matter content from the inverted elastic parameters according to their empirical relations. Zhao et al. (2015) and Verma et al. (2016) propose machine-learning approaches for organic-matter content prediction.

In this paper, we first analyze the geologic background and well logs of the target Lower Silurian shale formation. Then, a new anisotropic AVO inversion based on a modified PP-wave AVO equation is applied to estimate three subsurface attributes (A, acoustic impedance [AI]; B, shear modulus proportional to the anellipticity parameter; and C, approximate horizontal P-wave phase velocity). We show that the shale formation can be effectively discriminated from the surrounding limestone and sandstone by using the crossplot of attributes A and C; and the gas-bearing shale can be easily identified through the combination of attributes C and B. Finally, we apply a statistical rock-physics method to the inverted attributes for quantitative interpretation of the shale-gas reservoir. This method combines the Bayesian rule and Markov random fields to ensure the accuracy and continuity of the interpretation result. The final results show an encouraging comparison with the borehole measurements.

Geology background and well-log analysis

The study area is located in the southern Sichuan Basin, southwest of China (Figure 1a). Well logs in this survey reveal geologic systems including Ordovician, Silurian, Devonian, and Permian, from bottom upward, as listed in Figure 1b. As one of the major shale-gas reservoir formations in China, the Longmaxi Group (S1l) is a lower Silurian-age unit of marine sedimentary rock consisting of overmature organic dark shale. It is typically found deposited on the Wufeng Group (O3w), an Upper Ordovician shale overlying limestone formations, and it was placed below the Shiniulan Group (S1s), which is carbonate platform facies consisting of limestone and calcareous sandstone. In this area, the S1l has a thickness of 300 m, and the O3w has a thickness of 5 m.

The logs of well 3 in this study area are analyzed in this section (Figure 2a). The well logs are calibrated with the inline seismic section at approximately CDP 400 (Figure 2b). The mineralogy log shows that S1l shale (2080–2390 m) has moderate clay contents and varied volume of quartz and carbonate minerals. The lower S1l formation (2210–2390 m) is organic-rich and gas saturated. It consists of large volume of quartz and has much higher porosity and organic matter content than the upper S1l formation (2080–2210 m). The rock properties of O3w are very similar to that of the lower S1l formation. Logs of C33 and C44 are calculated using the logs of velocities and density. The log of the anisotropy parameter γ is calculated using C44 and C66 measured in the borehole. The anisotropy parameters ϵ and δ are not measured directly in the borehole, so we apply a rock-physics-based method to predict the anisotropy parameters for the shale formation (Zhang et al., 2017; Zhang, 2019). Logs of mineralogy, porosities of clay-bound water, free water, and free gas, and organic matter content are used to model all five effective elastic constants of the VTI shale by using a combination of varied anisotropic effective-medium theories, e.g., Backus averaging, self-consistent approximation theory, and the differential effective medium theory. The predicted anisotropy parameters are demonstrated by the agreement between the predicted γ and the measured γ in Zhang (2019). The shale formation has generally lower elastic moduli (lower velocities and densities) than upper and lower limestone formations, but it shows much higher magnitude of anisotropy than surrounding formations. Zhang (2019) analyzes major factors of intrinsic anisotropy of the Lower Silurian shale: (1) The magnitude of anisotropy is mainly determined by clay contents and also increases with the clay-bound water porosity, and (2) the presence of organic matter and pores of free water/free gas has a negative contribution to anisotropy of this overmature shale because the former are observed to distribute randomly within the rock matrix, and the latter have varied geometry and a wide range of orientations.

Simultaneous AVO inversion of three attributes for the VTI media

Simultaneous AVO inversion method

In this section, we show the influence of anisotropy in seismic amplitude and wepresent a modified reflection coefficient approximation of VTI media. The elastic stiffness tensor of a VTI medium is  
Thomsen (1986) defines anisotropy parameters as  
Therefore, elastic stiffness coefficients can be expressed by using vertical velocities, bulk density, and anisotropy parameters as  
where VP0, VS0, and ρ are the P- and S-wave vertical velocities and the bulk density, respectively, and ε, γ, and δ are the P-wave, S-wave, and anelliptical anisotropy parameters, respectively, proposed by Thomsen (1986).
Considering a planar interface separating two VTI layers, a linear approximation of the PP-wave reflection coefficient can be expressed as (Zhang et al., 2018, 2019)  
where θ is the incidence phase angle, and it can be approximated to the average angle θ of the incidence phase angle and the transmission phase angle under the assumption of weak impedance contrast and small angle; A=ρVP0, B=ρVS02eσ4, and C=VP0eε are the three attributes to be inverted for, σ=(VP0/VS0)2(εδ) is the effective parameter introduced by Tsvankin and Thomsen (1994), K=(2VS0/VP0)2, and Δ stands for the difference value for two layers across the interface. Attribute A represents AI, the attribute B is a shear modulus weighted by eσ4 in terms of an anellipticity parameter (εδ), and attribute C=VP0eε can be equivalent to the horizontal P-wave velocity. This approximation is compared with the exact solutions of anisotropic media (Schoenberg and Protázio, 1992) and isotropic media (Aki and Richards, 1980), respectively, as shown in Figure 3, in which the isotropic linear approximation (Aki and Richards, 1980) is also plotted in.

Table 1 shows the elastic parameters of four single-interface models belonging to different AVO classes. All of the model parameters are derived from well logs in the study area. The class I and class II models represent typical shale/sand interfaces; and class III and class IV represent the interface of the top of the gas-bearing shale (2360 m in well 3) and the top of the shale formation (2080 m in well 3), respectively. The anisotropy parameters of the shale formation are ε=0.15, δ=0.05, and γ=0.2 according to the well logs. The exact isotropic equation deviates from the exact anisotropic equation and such deviation increases with angle. The deviation for classes I, II, and III is larger than that of class IV. The anisotropic approximation (equation 4) has better accuracy than the exact isotropic equation for classes I, II, and III, and it has similar accuracy to the exact isotropic equation for class IV. The differences between approximations to their corresponding exact equations for class I and class IV are larger than those of class II and class III, and such differences are mainly caused by the strong impedance contrast across the interface. There is a huge difference between the isotropic approximation and the exact isotropic equation for class I. This error is also caused by the strong impedance contrast across the interface, but it happens to make the isotropic approximation have a more similar value to the exact anisotropic equation than the other three equations within 0°–40°. Therefore, simultaneous inversion based on the new anisotropic equation is more appropriate than the conventional isotropic inversion. The anisotropic inversion of the three attributes (A=ρVP0, B=ρVS02eσ4, and C=VP0eε) is no more difficult than the conventional isotropic inversion, whose ill posedness can be effectively reduced by an appropriate regularization. In the next section, we shall discuss the inversion scheme and demonstrate it using real seismic data.

Equation 4 consist of model parameters in the differential form of the natural logarithm. Therefore their original values, instead of reflectivities, can be directly estimated by using an inversion scheme similar to Stolt and Weglein (1985) and Buland and Omre (2003). The model parameters can be solved in a Bayesian framework by minimizing a quadratic misfit function J regularized by an a priori model of Gaussian distribution (Kennett et al., 1988; Wang and Houseman, 1994; Simmons and Backus, 1996; Tarantola, 2005)  
where g is the nonlinear forward-modeling operator, d is the data vector consisting of seismic amplitude as a function of incident phase angle at each time sample, m is the model vector consisting of the three model parameters to be solved, mprior is the prior model parameter vector, Cd is the data covariance matrix and could be practically calculated by assuming uncorrelated noise as Cd=σn2I, σn is the noise standard deviation, Cm=[CAACABCACCBACBBCBCCCACCBCCC] is the model covariance matrix that can be obtained from real logs or rock-physics relations, and Ckl (k and l refer to attributes A, B, and C) is the covariance of inverted parameters. Given an initial model mprior, the forward modeling can be approximated as  
where G is the linear mapping operator from model to data. It includes the effect of incident angles and wavelets, and it is given by (Downton and Ursenbach, 2006)  
where W is the convolution matrix containing wavelets of different angles, F is the sensitivity matrix of Fréchet derivatives with respect to the model parameters, and D is the differential operator of the model parameters in adjacent layers (Stolt and Weglein, 1985). The model parameter can be recovered iteratively as  
where Δd=dg(mprior) is the data residual and μ proportional to σn is the trade-off parameter of the regularization term.

Synthetic tests

Well logs (VP0, VS0, density, Thomsen’s ϵ and δ ) transformed from the depth to the time domain are used to generate the synthetic gathers to test the inversion scheme. Each gather consists of five traces from 10° to 50°, and the reflection coefficients are calculated using the exact solution in Schoenberg and Protázio (1992). Figure 4 shows the synthetic angle gathers with different signal-to-noise ratios (S/Ns). The noise variance σn in synthetic tests is 0.0003, 0.001, and 0.005 for the case of S/N = 2, 1, and 0.5, respectively, and μ can be estimated accordingly. The real logs and inverted attributes are compared as shown in Figure 5. The real logs are smoothed by a moving average filter with a 40-point span to generate the initial models. The inversion stability is analyzed using the relative error (RE) and crosscorrelation coefficient (CC) between the true model and the inversion result (Table 2). RE is calculated as i=1N|(minv(i)mreal(i))/mreal(i)|, where mreal(i) and minv(i) refer to the true and inverted parameters of the ith time sample, respectively. All three attributes (A=ρVP0, B=ρVS02eσ4, and C=VP0eε) can be satisfactorily inverted from the synthetic data with high S/N (S/N>2). REs increase with the decrease of S/N, whereas CCs decrease with the S/N. All three inverted attributes tend to be very patchy when using noisy data (S/N<1), and those of attribute B have more error and larger uncertainty than those of attributes A and C. The above analysis implies that the inversion scheme can be feasible for a field-data application using an S/N seismic data set with S/N >1.

Field-data inversion

In the field-data application, a series of processing procedures are used to enhance the S/N and to preserve the true reflection amplitudes, which is important for AVO inversions. Multiples and ground roll of original seismic gathers are first attenuated to enhance the S/N and preserve amplitudes of near-offset data. Double-parameter scanning, nonhyperbolic moveout, and anisotropic Kirchhoff time migration are used to flatten and image the reflection amplitude at far offset. The NMO stretch is also minimized by using the wavelet-stretch correction. A practical inversion workflow is performed on the common-image-point gathers: (1) transformation of the prestack seismic gathers from the offset domain to the angle domain, (2) seismic data within certain angle intervals are stacked to construct the constant-angle sections, (3) estimation of the mixed-phase wavelets for each of the constant-angle sections based on a high-order statistic method (Lü and Wang, 2007), (4) building initial models using the smoothed well logs (calibrated in the time domain) and picked horizons, and (5) iterative linear inversion for the three attributes (A, B, and C), and the calibrated well logs are used to calculate the model covariance for regularization. Steps 1 and 4 are performed using the Hampson-Russell software. The constant angle sections and estimated mixed-phase wavelets are shown in Figures 6 and 7. The inversion results are shown in Figure 8. There is a good agreement between the inversion results and the corresponding logs (Figure 9a9c), and the synthetic angle gather generated using the inverted parameters and real seismic angle gather are also compared in Figure 9d and 9e. There are several differences between the inversion results and log data within 1280–1380 ms. These differences correspond to the differences between the seismic data and calibrated well logs. Few fluctuations are found in the log data, and, thus, no reflection events are found in the synthetic data gather (Figure 4a) at the same time interval. The shale formation has lower impedance than the upper and lower formations and thus can be roughly interpreted within the time window from 1260 to 1460 ms at the well location (CDP 400). The shale-gas reservoirs are approximately within the range of 1440–1460 ms, in which the inverted attributes show a lower value compared with the upper shale layers. In the next section, we apply a statistical-rock-physics-based method to the inversion results for quantitative interpretation of this shale-gas reservoir.

Shale-gas reservoir characterization using statistical-rock-physics-based quantitative interpretation

The interpretation capability of three inverted attributes in shale-gas characterization is discussed using log data. The crossplot of A-C, representing the AI and VP0eε, respectively, is compared with that of AI-VP0 in Figure 10. It is easy to identify the shales samples (clay volume >15%, AI<14km/s*g/cm3, C<6.0km/s) from surrounding carbonates in the crossplot of A and C. This implies that attributes A and C are useful in shale formation prediction because shales show a distinct A-C relation from upper and lower formations due to its generally pronounced VTI properties. It is difficult to discriminate the shale from the surrounding formations by using the crossplot of AI and VP0 because there is a good linear relation between them. Once the shale formation is identified, we need to identify the organic-rich shale sections within the whole shale formation. The crossplots of C-B and the VP0-shear modulus (ρVS02) are compared in Figure 11. Only the logs within the shale formations (2080–2390 m) are used for this analysis. The gas-bearing organic-rich shale is identified as those of a high volume of kerogen and high porosity. It can be seen that this target is more easily identified by using the crossplot of C-B than using VP0-ρVS02, which is useful in conventional oil and gas reservoir interpretation. Besides, in fact, interpretation of field data using isotropic inversion results (e.g., ρVS02 and VP0) may be less accurate than that of the anisotropic inversion results (attributes B and C) because the former can be biased by the anisotropy. Therefore, the shale-gas reservoir can be characterized by a two-step interpretation using the inverted attributes: The shale formation is first discriminated from surrounding limestones by using the crossplot of the AI and the anisotropic P-wave velocity; then the gas-bearing shale is identified through the combination of the anisotropic P-wave velocity and the anisotropic shear modulus and the AI.

Mapping the organic-rich gas-bearing shale distribution from seismic inversion results is of great importance in shale-gas reservoir characterization and further development. By using crossplots of well-log data, we have shown that the inverted attributes can be useful for shale-gas reservoir characterization. In this section, we apply a quantitative interpretation method proposed by Wang et al. (2019) to the inverted attributes as shown in Figure 12. The quantitative interpretation workflow uses statistical rock-physics methods as discussed in Mukerji et al. (2001) and Avseth et al. (2005), and the Markov random field is combined with the Bayesian rule to improve the continuity and accuracy of the interpreted results. First, shale, organic-rich (gas-bearing) shale and limestone are defined as three lithofacies to be classified. Well logs of clay content, kerogen content, and porosity are used as training data to classify the three lithofacies based on the crossplots of A-C and C-B. Second, the scale of well logs is expanded to that of inverted seismic results with empirical value λ/8 (λ represents the wavelength of seismic data) by using the Backus average (Backus, 1962), and correlated Monte Carlo (CMC) simulations (Avseth et al., 2005) are used to extend number of data points for each lithofacies. Figure 13 shows the comparison between CMC result and raw well-log data of attributes A and C corresponding to limestone and shale. Figure 14 shows the comparison between CMC result and raw well-log data of attributes C and B corresponding to shale with low organic matter content and organic-rich shale. Third, the kernel-based, nonparametric, probability-density estimation is performed to construct the 2D conditional probability density functions (CPDFs) of attribute pairs (A-C and C-B) for varied lithofacies (Figure 15a and 15b). The elements of a Bayesian confusion matrix represent the conditional probability of being the true lithofacies given a predicted lithofacies are shown in Figure 15c and 15d. The conditional probabilities in both confusion matrices have very high values. Fourth, the Bayesian rule and Markov random fields are introduced in corresponding inverted seismic attributes for lithofacies classification. In this step, an initial lithofacies classification is estimated according to the CPDFs of training data from well logs. The most probable lithofacies classification is obtained iteratively according to the Bayesian rule, and a Markov random field is applied in the estimation of the prior probability in each iteration to maintain high spatial and vertical continuity. Figures 16 and 17 show the interpretation results of limestone, shale, and organic gas-bearing shale formations from inverted attributes. The prediction results in the well location agree well with the lithofacies classification extracted from the well log.


This paper focuses on the shale-gas reservoir characterization by using an AVO inversion method for VTI media and a statistical-rock-physics-based quantitative interpretation method. The AVO inversion for VTI media is based on a modified anisotropic PP-wave reflection coefficient consisting of three attributes. Three attributes included in the new equation have explicit physical and can be well-inverted from prestack seismic data using the simultaneous inversion in a Bayesian framework. Practical procedures of conventional isotropic inversion, such as wavelet estimation, well tie, and initial-model building, can be used in this inversion. Note that independent anisotropy parameters (ϵ and δ) are not directly inverted by using this inversion method. They could be recovered from inverted attribute B and attribute C through a stepwise inversion scheme, and this may be a potential use of this study. Inverted attributes show good interpretation capability in shale-gas reservoir characterization: Shales can be easily classified according to the nonlinear relation in the crossplot of attributes A and C, and the organic-rich gas-bearing shale is identified through the combination of attributes C and B. A statistical rock-physics workflow is applied to the inverted attributes for quantitative seismic interpretation. It combines the Bayesian rule and Markov random fields to ensure the accuracy and continuity of the interpretation result. Comparison between the prediction results and well logs demonstrates the feasibility of the inversion and quantitative interpretation approaches.


This study is sponsored by National Science and Technology Major Project (grant no. 2017ZX05018005) and CNPC science research and technology development project (grant no. 2019A-3308).

Data and materials availability

Data associated with this research are available and can be obtained by contacting the corresponding author.

Feng Zhang received a bachelor’s degree in electrical and mechanical engineering from Beijing Institute of Technology and a master’s degree in signal processing and communications from Edinburgh University. He started his Ph.D. research in reservoir geophysics at Imperial College London in 2006, and he received a doctorate degree from Imperial College in 2010. He is a professor at the China University of Petroleum (Beijing). He attended the Edinburgh Anisotropy Project as a post-doctoral research associate. In 2011, he joined the China University of Petroleum (Beijing). He is an associate editor of the Journal of Geophysics and Engineering. His research interests include seismic anisotropy, seismic inversion, rock physics, and reservoir characterization.

Lin Wang received a bachelor’s degree in applied geophysics from the China University of Petroleum (East China) and a master’s degree in geophysics from the China University of Petroleum (Beijing). He has been a Ph.D. candidate at the China University of Petroleum (Beijing) since 2016. His current research interests include statistical rock physics, shale gas reservoir characterization, and seismic anisotropy.

Xiang-yang Li received a bachelor’s degree in geophysics from Changchun Geological Institute, a master’s degree in applied geophysics from East China Petroleum Institute, and a doctorate degree (1992) from The University of Edinburgh. He is a professor at the China University of Petroleum (Beijing) and the project leader of Edinburgh Anisotropy Project. He has been working for British Geological Survey since 1991. He joined the China University of Petroleum as a professor in 2008. His current research activity includes multicomponent seismology, seismic data processing and interpretation, and seismic anisotropy.

Freely available online through the SEG open-access option.