Hydrocarbon Identification and Bedding Fracture Detection in Shale Gas Reservoirs Based on a Novel Seismic Dispersion Attribute Inversion Method

Prediction of hydrocarbon enrichment and natural fractures is significant for sweet spot characterization in shale gas reservoirs. However, it is difficult to estimate reservoir properties using conventional seismic techniques based on elastic and isotropic assumptions. Considering that the viscoelastic anisotropic model better represents organic shale, we propose a new seismic inversion method to improve shale gas characterization by incorporating the anisotropic reflectivity theory in the frequency-dependent inversion scheme. The computed P-wave velocity dispersion attribute DP evaluates the hydrocarbon enrichment by estimating the inelastic properties of shale associated with organic materials. The inverted anisotropic dispersion attribute Dε detects the development intensity of bedding fractures using frequency-dependent anisotropy owing to wave-induced fluid flow in parallel fractures. Synthetic tests indicate that DP can robustly estimate shale attenuation and Dε is sensitive to the frequency-dependent anisotropy of shale. The results are validated by reservoir properties measured in gas-producing boreholes and rock physical modeling analysis, supporting the applicability of the dispersion attributes for hydrocarbon identification and bedding fracture detection. The predicted hydrocarbon enrichment and the development of bedding fractures correlate with the structural characteristics of the shale formation. The depth-related shale properties can be described by improving the geological understanding of the study area. Finally, favorable areas with high hydrocarbon enrichment and extensive development of bedding fractures are identified by simultaneously considering high DP and Dε anomalies, providing essential information for predicting potential shale gas reservoirs. A novel seismic inversion method for anisotropy dispersion attributes is proposed P-wave velocity dispersion attribute is used to identify hydrocarbon enrichment in shale Anisotropic dispersion attribute is used to detect bedding fractures in shale A novel seismic inversion method for anisotropy dispersion attributes is proposed P-wave velocity dispersion attribute is used to identify hydrocarbon enrichment in shale Anisotropic dispersion attribute is used to detect bedding fractures in shale


Introduction
Hydrocarbon enrichment and the development of natural fractures are essential for identifying sweet spots in shale gas plays. Total organic carbon (TOC) has been used to determine gas generation in the organic shale of the Longmaxi Formation (Zou et al. 2015;He et al. 2017;Zhang et al. 2020a). Horizontal bedding fractures, accounting for most natural fractures in the Longmaxi shale, provide storage spaces and migration channels for gas enrichment (Zhang et al. 2020b;Xu et al. 2021) and are favorable for hydraulic fracturing (Zhou et al. 2020). Using seismic methods to identify hydrocarbon enrichment and predict bedding fractures can provide valuable information for locating promising shale gas reservoirs.
As density is considered an indicator of TOC in the Longmaxi shale (Li and Yin 2015), current seismic methods for identifying the variations in TOC and gas content mainly estimate density using improved pre-stack elastic inversion . However, the sensitivity of seismic reflection to rock density is weak, limiting robust density estimates using pre-stack seismic data. Moreover, the positive correlation between TOC and density may yield good results in some areas but cannot be applied universally, limiting its application for TOC estimates. On the other hand, current research is focused on predicting vertical fractures in shale using wide-azimuth seismic inversion as they exhibit particular azimuthal seismic signatures that can be conveniently extracted (Bachrach 2015;Kumar et al. 2018;Luo et al. 2019). However, because of intense tectonic uplift and overpressure formed by gas generation, massive bedding fractures are observed in the Longmaxi shale that can accumulate and preserve natural gas (Zhang et al. 2020b;Xu et al. 2021). Thus, it is necessary to develop methods to predict bedding fractures in shales. However, seismic methods are rarely used for such predictions because seismic signatures of horizontal fractures are not as distinct as those of vertical fractures.
A rock physics model that considers viscoelasticity and anisotropy is more appropriate to describe the characteristics of shale (Carcione 2001;Carcione et al. 2011). Viscoelasticity describes the inelastic properties of organic shale associated with the organic matter composed of kerogen and hydrocarbon-saturated pores. Several factors can cause vertical transverse isotropy (VTI) in shale, including the anisotropy of clays, layered distribution of organic matter, laminated fabric of the solid matrix, and bedding fractures generated due to overpressure or developed along weak planes in the shale matrix. Several studies have been conducted on the effect of the anisotropic clay minerals and laminated organic material and fabric on shale anisotropy (Vernik and Nur 1992;Vernik and Landis 1996;Carcione 2000Carcione , 2001Carcione et al. 2011;Sayers 2013;Sayers et al. 2015). Bedding fractures developed in the shale matrix are another important factor that controls shale VTI anisotropy. Vernik (1994) and Vernik and Liu (1997) experimentally investigated the stress-related anisotropy of organic shale, which could be attributed to the development of bedding fractures related to overpressure during maturity of the organic matter. The closing and opening of aligned fractures at different effective stress fields can cause variations in shale anisotropy. Recent experimental studies have shown that the Longmaxi shale exhibits a substantially variable anisotropy with loading stress fields (Liu et al. 2019b), indicating the presence of bedding fractures, which has been further supported by the observations from cores and outcrops (Zhang et al. 2020b;Xu et al. 2021).
Based on the geological "sweet window" theory , it has been suggested that owing to intense tectonic activity in the Longmaxi Formation, bedding fractures exhibit a more heterogeneous spatial distribution than the lamination fabric of shale.
According to the equivalent media theories of penny-shaped cracks (Hudson 1981(Hudson , 1988 and linear-slip fracture models (Schoenberg 1980;Schoenberg and Sayers 1995), a few aligned fractures saturated with hydrocarbon fluids would lead to intense transverse isotropy (TI) because of the remarkably contrasting compliance between soft fracture inclusions and a solid matrix composed of minerals. Therefore, it is expected that heterogeneity of the shale VTI anisotropy in the study area could be primarily because of the variably intense bedding fracture development, with minor inputs from other factors, such as clay minerals or organic matter. Thus, anisotropic parameter estimates using seismic methods can be used to detect bedding fractures. Chapman (2003) suggested that seismic waves propagating in media with parallel fracture systems could experience anisotropic dispersion and attenuation owing to waveinduced fluid flow, leading to the frequency-dependent anisotropy of the rock. Furthermore, it can reasonably be assumed that clay minerals and shale laminae fabrics are elastic and do not produce anisotropic dispersion. Additionally, laboratory experimentations have indicated that the organic matter in the Longmaxi shale is usually distributed in scatter patches (Li et al. 2018), negligibly influencing shale VTI anisotropy. Therefore, it is feasible to detect bedding fractures using anisotropic seismic dispersion. Wilson et al. (2009) used the frequency-dependent amplitude variations versus offset (FD-AVO) scheme to estimate the seismic dispersion attributes. Subsequently, FD-AVO methods derived from various linearized approximations of PP-wave reflection coefficients were proposed and applied to shale gas reservoirs (Pang et al. 2018;Liu et al. 2019a). However, the main limitation of these methods is assuming the medium to be isotropic, which is inappropriate for anisotropic shale. A coupled issue is that most current FD-AVO techniques only estimate P-wave velocity dispersion, neglecting valuable information at far offsets in recorded seismic reflections, which can be used to estimate anisotropic dispersion for bedding fracture detection.
In this study, we propose a novel seismic inversion method to compute P-wave velocity and anisotropic dispersion attributes by extending the conventional frequency-dependent inversion scheme using anisotropic reflectivity theory. The computed dispersion attributes are used for hydrocarbon identification and bedding fracture detection in shale gas reservoirs. First, the proposed inversion method is presented in detail. We then discuss the frequency-dependent reflection behaviors related to the inelastic and anisotropic properties of shale. Subsequently, the sensitivity and robustness of the proposed dispersion attributes to inelastic properties of shale are tested using synthetic data. Finally, the developed method is applied to field data to evaluate the hydrocarbon enrichment and development degree of bedding fractures in the Longmaxi shale, and favorable shale gas zones are identified based on the computed results.

Seismic Anisotropic Dispersion Attribute Inversion
Schematic diagrams of the seismic reflection model are shown in Fig. 1, where a viscoelastic VTI shale is surrounded by two elastic isotropic media, representing the strata of the Longmaxi shale gas reservoir. The shale is overlain by sandstone and underlain by limestone. The plot on the right of Fig. 1 shows the relative variation in the P-wave impedance of the three layers. Rüger (1997) derived the PP-reflection coefficient for an interface separating two elastic VTI layers as follows: where V P and V S represent vertical P-wave and S-wave velocities, respectively. Z and G are the vertical P-wave impedance and vertical S-wave modulus. ε and δ are the Thomsen anisotropy parameters. Symbols Δ and − indicate the difference and average of the corresponding properties across the boundary, respectively.
Based on the following relationship Equation (1) can be rearranged as

By setting
Equation (3) can be expressed as (1) Fig. 1 Schematic diagrams of the model where a viscoelastic VTI shale is surrounded by two elastic isotropic media. Relative variation of the P-wave impedance for the three media is displayed, representing the model for the Longmaxi shale formation According to the FD-AVO framework, we extended Eq. (5) to the frequency domain where the elastic properties and anisotropy parameters are assumed to be frequency dependent owing to the viscoelastic anisotropy of shale. The density can be considered independent of the frequency. Frequency-dependent rock properties cause reflection coefficients to vary with frequency, which is the core of the FD-AVO scheme.
Using the first-order Taylor series, we expanded Eq. (6) and defined the PP-reflection coefficient at the reference frequency f 0 as We obtained the difference in the PP-reflection coefficient with respect to that at the reference frequency f 0 For field data applications, we implemented spectral decomposition to Eq. (8) and incorporated the effect of the wavelet with the spectrum W( f) where ΔS(θ, f) denotes the decomposed spectra of the pre-stack angle gather. D P , D G , D δ , and D ε are the dispersion attributes related to P-wave velocity, shear modulus, and anisotropy parameters δ and ε, respectively.
For the decomposed spectra ΔS(θ, f) of a pre-stack angle gather with n incidence angles and m frequencies, Eq. (9) can be expressed as a matrix to invert the dispersion attributes in Eq. (10) Equation (11) can be expressed in a simplified form as: where d represents the decomposed spectra ΔS(θ, f) and G is the coefficient matrix in Eq. (11). Then, the dispersion attributes in Eq. (12) can be computed using the least square method where I denotes the identity matrix and σ represents the damping factor. Cohen (1995) proposed the smooth pseudo-Wigner-Ville distribution method to improve the traditional Wigner-Ville distribution for the decomposition of nonstationary seismic signals as follows:

Method for Spectral Decomposition Estimate
where X(t) represents the analytical signal of the seismic response x(t) based on the Hilbert transform. S(t, f) represents the decomposed spectrum of X(t). Parameters v and τ are the time delay and frequency offset, respectively. g and h represent the window functions in the frequency and time domains, respectively, and are proposed to reduce the cross-term interfaces. We employed the smooth pseudo-Wigner-Ville distribution method to obtain the decomposed spectra of seismic data in the inversion using Eq. (11). (11)

Viscoelastic Anisotropic Model of Organic Shale
Organic shale can exhibit inelastic properties owing to the intrinsic attenuation of the organic mixture composed of kerogen and hydrocarbon-saturated pores (Carcione 2000(Carcione , 2001. Seismic anisotropy of organic shale has diverse geological origins, including the preferred orientation of clay minerals and organic matter, laminated fabric, and bedding fractures (Vernik and Nur 1992;Vernik and Liu 1997;Carcione et al. 2011;Sayers 2013;Sayers and den Boer 2018). Moreover, parallel fractures can cause anisotropic dispersion and attenuation owing to wave-induced fluid flow (Chapman 2003). Therefore, it is assumed that organic shale simultaneously exhibits inelastic and anisotropic properties. However, complex attenuation and anisotropy mechanisms limit the development of sophisticated rock physics models for describing the relevant characteristics of organic shale. In this sense, the viscoelastic anisotropic medium theory (Carcione 1997) can provide an equivalent representation of organic shale, where frequency-dependent complex stiffnesses representing VTI anisotropy are expressed as follows: where c* 11, c* 33, c* 55, and c* 13 are the stiffnesses at high-frequency elastic limit. K and D have the following forms: The Zener complex moduli in Eq. (15) are defined as where v = P and Q represent P-and S-wave relaxation mechanisms, respectively.
In Eq. (17), the relaxation time of stress and strain are calculated as follows: where Q P and Q S are the quality factors describing the attenuation of P-and S-wave, respectively. τ 0 is the relaxation time associated with rock properties. The inverse of τ 0 defines the characteristic frequency f 0 corresponding to the attenuation peak. In this study, similar to the study by Carcione (2001), all modeling examples considered a characteristic frequency f 0 = 35 Hz, simulating a situation where velocity dispersion and attenuation occur at the seismic frequency.
The frequency-dependent Thomsen anisotropy parameters for VTI shale can be derived from the stiffnesses in Eq. (15) as follows:

Seismic Modeling by Integrating Viscoelastic Anisotropic Theory and Propagator Matrix Method
For P-wave incidence on a viscoelastic VTI medium surrounded by two half-spaces, Carcione (2001)  where i P represents a vector related to the incident P-wave. A 1 and A 2 are matrices associated with the properties of the upper and lower half-spaces, respectively. B denotes the matrix of the embedded VTI layer. PMM provides a practical method that incorporates inelastic mechanisms in seismic response simulation because the rock physical model and PMM are both conducted in the frequency domain (Guo et al. 2015a, b). Thus, the dynamics information of the seismic reflections can be reasonably simulated. After calculating the frequency-dependent reflection coefficient (R PP ) using Eq. (20), the spectrum of the reflected waves (U f ) can be computed when the spectrum of the incident wavelet (W f ) is known.
The seismograms (u t ) can be computed from the reflection spectrum (U f ) using the inverse Fourier transform

Velocity Dispersion and Attenuation of Viscoelastic Anisotropic Shale
In the schematic diagrams of the seismic reflection model (Fig. 1), vertical P-and S-wave velocities of the Longmaxi shale at high-frequency elastic limit were set to V P = 4117 m/s and V S = 2300 m/s, and the density was ρ = 2455 kg/m 3 according to logging data. The properties of the upper sandstone were V P = 4250 m/s, V S = 2360 m/s, and ρ = 2640 kg/m 3 , and those of the lower limestone were V P = 5849 m/s, V S = 3128 m/s, and ρ = 2721 kg/m 3 . The thickness of the shale gas reservoir was set at 40 m.
We analyzed the frequency-dependent shale properties for different anisotropic parameters and quality factors based on the effective medium theory described in Sect. 2.3.1. The first set considered three shale models with different anisotropic parameters and constant quality factors. Anisotropic parameters (ε, γ, δ) vary from (0.25, 0.30, 0.22), (0.15, 0.18, 0.12) to (0.05, 0.06, 0.02), corresponding to the decreasing degree of anisotropy. Relationships between the anisotropic parameters in the three models were determined from core measurements of the Longmaxi shale in the laboratory (Deng et al. 2014. ε exhibits a positive linear correlation with γ, and γ is higher than ε. δ was set according to the measured ε. The quality factors were kept constant at Q P = 20 and Q S = 15 for all three models, corresponding to a relatively strong attenuation. Numerous combinations of stiffnesses c* 11 and c* 33 can lead to the same anisotropic parameter ε. To analyze the effect of anisotropic parameters on frequency-dependent properties of shale, we maintained a constant vertical stiffness c* 33 and assumed that the change in ε is attributed to the horizontal stiffness c* 11, thereby avoiding the possible confusion of varying c* 33 owing to different quality factors in the following models. Figure 2a, b, and c shows the calculated frequency-dependent c 11, c 33, and P-wave anisotropic parameter ε, respectively. At the high-frequency elastic limit, ε equals the values set for the three models. The second set of models investigates the effect of quality factors, maintaining invariant VTI anisotropic parameters (ε = 0.15, γ = 0.18, δ = 0.12). We varied the P-and S-wave quality factors (Q P , Q S ) from (20, 15), (50, 30), to (100, 80), corresponding to a decreasing Fig. 2 Frequency-dependent stiffnesses of a c 11 , b c 33 , and c P-wave anisotropic parameter ε, for three shale models with variable high-frequency elastic anisotropic parameters. In the three cases, vertical P-and S-wave velocities at high frequency are set to 4117 and 2300 m/s. Quality factors are kept constant as Q P = 20 and Q S = 15, respectively attenuation in shale. Increasing attenuation (decreasing quality factors) reduces the c 11 and c 33 values at low frequencies, enhancing their degree of frequency-dependence of c 11 and c 33 ( Fig. 3a and b). The impact of attenuation on parameter ε is subtle (Fig. 3c).

PP-Reflection Coefficients for Viscoelastic Anisotropic Shale Layer
We investigated the effect of anisotropic parameters and quality factors on the frequencydependent reflection coefficients R PP using the PMM for the Longmaxi shale model shown in Fig. 1.
The effect of the anisotropic parameters on R PP , corresponding to the three shale models in Fig. 2, is illustrated in Fig. 4. Plots in the left column correspond to R PP from the top interface of the shale layer, and those in the right column correspond to the bottom interface. As indicated by the relative variation in P-wave impedance (Fig. 1), the top interface of the shale layer presents a type IV AVO response, and the bottom exhibits a type I AVO response. R PP at different frequencies were separated from each other for both the top and bottom reflections because of the frequency-dependent inelastic properties of shale. The separation of the AVO responses in Eq. (6) at different frequencies is fundamental for the feasibility of the FD-AVO method. The effect of shale anisotropy on R PP became evident Fig. 3 Frequency-dependent stiffnesses of a c 11 , b c 33 , and c P-wave anisotropic parameter ε, for three shale models with variable quality factors Q P and Q S . In the three cases, vertical P-and S-wave velocities at high frequency are set to 4117 and 2300 m/s. High-frequency elastic anisotropic parameters ε, γ, and δ are kept constant as 0.15, 0.18, and 0.12, respectively at incidence angles > 20°, suggesting that the variations in anisotropy were detectable at far offsets.
The impact of the quality factors on the frequency-dependent R PP is shown in Fig. 5. The three cases correspond to shale models with variable quality factors and constant anisotropy (Fig. 3). R PP curves were more separated for low-quality factors (high attenuation), which forms the base for conducting dispersion attribute inversion to estimate the inelastic properties of shale.

Sensitivity Analysis of Dispersion Attributes to Shale Properties
Sensitivity of the dispersion attributes D P and D ε to the inelastic and anisotropic properties of shale was investigated using synthetic seismic data of the Longmaxi shale model (Fig. 1), which were computed using the PMM method (Eqs. 20-22). Figure 6a, b, and c shows the computed AVO responses for the three shale models with different anisotropic parameters and constant quality factors Q P = 20 and Q S = 15 (Fig. 2). The incident Ricker wavelet had a dominant frequency of 35 Hz. Shadowed windows denote the reflections corresponding to the target shale layer. The modeling results indicate that the impact of shale anisotropy on seismic reflections becomes evident for incidence angles > 20°, which is consistent with the reflection coefficient analyses (Fig. 4). Figure 6d and e shows D P and D ε computed using the inversion scheme proposed in our study. In Fig. 6d, the inverted D P shows no variations, agreeing with the shale models that have constant quality factors. D ε presents distinct changes for the three models with varied anisotropy parameters (Fig. 6e). The results indicate that D ε is sensitive to the variations in the frequency-dependent VTI anisotropy of shale. Figure 7 shows the corresponding results for the three shale models with different quality factors and constant anisotropy. The effect of shale inelasticity is evident on the seismic reflections at all incidence angles (Fig. 7a, b, and c) that exhibit variations in the amplitude and phase of the reflected waves. D P represents a higher anomaly for shale with higher attenuation (Fig. 7d) and can be used to evaluate the inelastic properties of shale. Simultaneously, the difference in the computed D ε is little (Fig. 7e) because the variation in the frequency-dependent anisotropic parameters is subtle for the three models (Fig. 3c).
The results in Figs. 6 and 7 justify the feasibility of using D P and D ε to estimate the inelastic and anisotropic properties of organic shale. D P is sensitive to variations in quality factors, and D ε can robustly reveal the changes in the frequency-dependent anisotropic parameters of shale. Therefore, the computed D P and D ε can be used for hydrocarbon identification and bedding fracture detection in shale gas reservoirs. D P evaluates hydrocarbon enrichment using P-wave velocity dispersion associated with organic matter and hydrocarbon-saturated pores. At the same time, D ε predicts bedding fractures by considering their frequency-dependent anisotropy associated with fluid saturation.

Datasets
We applied our proposed method for the characterization of the Longmaxi shale gas reservoir. Figure 8a shows the regional geology map, where the red rectangle indicates the study area. Intense local tectonic movements formed a series of NE-SW-trending anticline and syncline belts. The traveling time of seismic reflections from the Longmaxi shale illustrates the tectonic relief of the target formation (Fig. 8b). Well A is located near the structural high of the Pingqiao Anticline Belt, and Well B is deployed on the slope of the Baima Syncline Belt. Figure 9a, b, and c displays horizontal slices of inverted seismic elastic properties of P-wave impedance (I P ), velocity ratio (V P /V S ), and density (ρ) of the Longmaxi shale gas reservoir, respectively. Figure 10a shows the seismic profile across Well A and Well B. The bottom of the Longmaxi shale is indicated by a yellow line. Owing to severe intense tectonism in the area, several anticlines and synthetics added significant variations to the burial depth of the Longmaxi shale formation. Figure 10b demonstrates the flattened seismic profile of the Longmaxi shale. Figure 11a, b, and c illustrates cross-well profiles of I P , V P /V S ratio, and ρ of the target shale, respectively.
The continuous distribution of I P over a vast area in the Longmaxi shale suggests that I P may be not an efficient indicator for reservoir characterization. Previous research indicates that shale density is negatively correlated with TOC and can be regarded as an indicator of hydrocarbon accumulation (Schmoker 1979;Alfred and Vernik 2013;Lim et al. 2017). This correlation has been utilized to characterize the Longmaxi shale gas reservoir . Shale around the gas-producing Well B exhibits a low density (Fig. 9b), implying a potential area for gas enrichment. However, the gas-producing Well A shows a high shale density, implying that density may not be a robust indicator of hydrocarbon enrichment in the study area. Furthermore, the use of V P /V S ratio (Figs. 9c and 11c) as a hydrocarbon indicator for the promising shale gas reservoir is challenging because the poroelastic behavior of shale is different from that of sandstone, making Gassmann's fluid Fig. 7 Computed seismic AVO responses for three models a, b, and c, with different anisotropic parameters. d Inverted P-wave velocity dispersion attribute D P . e Anisotropic dispersion attribute D ε substitution theory inapplicable. Nevertheless, the V P /V S ratio is critical for evaluating the shale brittleness index in hydraulic fracturing (Grieser and Bray 2007;Guo et al. 2012).
The above analyses suggest that inverted seismic elastic properties may not serve as robust indicators for hydrocarbon identification. Furthermore, the assumption of an isotropic medium in the corresponding methods limits their applications to fracture prediction. Therefore, practical methods that consider the inelasticity and anisotropy of shale are necessary for effective hydrocarbon identification and fracture detection in shale gas reservoirs.

Seismic Dispersion Attributes and their Implications of Shale Reservoir Properties
The proposed inversion method for seismic dispersion attributes was applied to the Longmaxi shale gas reservoir in the study area shown in Fig. 8. Figure 12a and b illustrates the horizontal slices of the computed D P and D ε of the Longmaxi shale gas reservoir. First, a large region on the Fenglai Syncline Belt, located northwest of the study area, is characterized by high D P and D ε anomalies and can be considered as a potential area for future shale gas exploration and development. However, shale gas reservoirs in this area have not been explored because of their considerable burial depth. Owing to the relatively higher D P and D ε in this area, the dynamic display range may require adjustment for detailed characterization. Moreover, to interpret the high dispersion anomaly in this region, the poroelastic Fig. 9 a Horizontal slices of P-wave impedance I P . b Density ρ. c Velocity ratio V P /V S for the Longmaxi shale behaviors of deep shale gas reservoirs should be modeled using rock physics methods, which is beyond the scope of this study. Thus, we focus on other promising shale gas zones with relatively shallower burial depths. Apart from the deep Fenglai Syncline Belt, three favorable gas zones on the Jinping and Pingqiao Anticline Belts and the Baima Syncline Belt can be distinguished from the D P and D ε ( Fig. 12a and b). In comparison, the Shuanghekou Syncline Belt presents low D P and D ε , which may not be considered as a potential area. For the three promising zones, the correlation between structures and dispersion attributes may have particular geological implications, such as depth-related reservoir properties of shale gas plays , which deserve further investigation based on an improved geological understanding of the local area. The Jinping Anticline Belt simultaneously exhibited apparent D P and D ε anomalies. This supports that hydrocarbon enrichment indicated by a high D P anomaly may positively correlate with the intensity of development of bedding fractures generated by overpressure owing to gas generation in this region (Xu et al. 2021), as suggested by a high D ε anomaly. Therefore, the Jinping Anticline Belt is a promising area for future exploration. The Longmaxi shale in the Pingqiao Anticline Belt exhibits a relatively high D P anomaly and a much higher D ε anomaly. In comparison, the shale in the Baima Syncline Belt presents a very high D P anomaly and a lower D ε value. The structural characteristics of dispersion attributes of the shale in the Pingqiao Anticline Belt and Baima Syncline Belt  Fig. 12 will be further investigated on the profiles of D P and a D ε across Well A and Well B, as shown in Fig. 13.
A high D P anomaly is observed in Well A and Well B, which is more significant around Well B (Fig. 13a). From the implication of D P proposed in our study, we can infer that the shale in Well B has more hydrocarbon enrichment than that in Well A. Table 1 shows the average reservoir properties (TOC, porosity, total gas content, and layer thickness) of the Longmaxi shale estimated from the data of the two gas wells. Values of all relevant reservoir properties are higher in Well B than those in Well A. Higher TOC, porosity, and gas saturation in shale correspond to higher organic content and cause a more evident inelasticity revealed by high D P values. Additionally, the shale is thicker in Well B (Table 1), causing a much more significant D P anomaly as can be expected. These results justify the applicability of the dispersion attribute D P for hydrocarbon identification.
Meanwhile, Well A has a much higher D ε anomaly than that in Well B (Fig. 13b). The D ε values suggest that the shale in Well A should have developed more intense bedding fractures than that in Well B (Fig. 13b). The Longmaxi shale in Well A is located on the structural high of the Pingqiao Anticline at a shallower burial depth (Fig. 10a) than that in Well B, which allows the development of massive bedding fractures along the mechanically weak surface of the shale rock owing to overburden off-loading and stress release (Zhang et. 2019(Zhang et. , 2020bXu et al. 2021). As shown by the formation micro-scanner image (FMI) logs (Figs. 14 and 15), the target shales in Well A developed more bedding fractures than those in Well B.  Additionally, we conducted a rock physical analysis to evaluate the development degree of bedding fractures for the shale in Well A and Well B using the modeling method of Guo et al. (2022). The computed Hashin-Shtrikman bounds were illustrated on the bulk modulus versus porosity cross-plot (Fig. 16). In the modeling, kerogen was treated as an inclusion in the shale matrix and is a part of the porosity. The properties of the different components of the rock are provided in Table 2. The corresponding averaged volumetric fractions of the constituents were estimated using the logging data. The superimposed data points for shales were extracted from the logging data in the two wells. Difficulty in accurately obtaining elastic properties of kerogen and clay minerals and possible variations in kerogen distribution in shale matrix may cause uncertainty in calculating the Hashin-Shtrikman bounds. Nevertheless, the distribution of data points in Fig. 16 can reveal the difference in the degree of bedding fracture development.
The shale in Well A presented a relatively lower bulk modulus for a specific porosity and tended to be close to the lower bound (Fig. 16). With respect to the Hashin-Shtrikman bounds, a lower bulk modulus close to the lower bound for a given porosity indicates a more fractured or cracked pore space, corresponding to increasing fracture density. Considering bedding fractures as the primary fracture type in the Longmaxi shale (Zhang et al. 2019(Zhang et al. , 2020b, we can infer a higher intensity of bedding fracture development in Well A, agreeing with the higher D ε values for the shale in Well A (Fig. 13b). Therefore, these results validate that the proposed D ε can be used for bedding fracture detection.

Application of the Dispersion Attributes to Predict Favorable Shale Gas Zones
As discussed above, the dispersion attributes, D P and D ε , have specific implications for reservoir properties in terms of hydrocarbon enrichment and bedding fracture intensity in organic shale and can be used to predict favorable shale gas reservoirs. Figure 17 illustrates the results estimated from Fig. 12a and b, where the thresholds of normalized D P and D ε values were set at 0.5 and 0.65, respectively. Apart from the potential area in the Fenglai Syncline Belt, three favorable zones were identified. However, the threshold values of D P and D ε were determined empirically, as the quantitative relationships between dispersion attributes and reservoir properties have not been investigated. The threshold should be optimized in future studies with sufficient geological and geophysical information.

Conclusions
The viscoelastic anisotropic model is considered to better represent organic shale. Accordingly, a novel seismic inversion method was proposed to compute the dispersion attributes to estimate hydrocarbon enrichment and bedding fracture development in shale gas reservoirs. The method incorporates the anisotropic reflectivity theory in the frequencydependent inversion scheme, which was proven to be more appropriate for shale gas characterization than conventional methods based on elastic and isotropic assumptions. Synthetic tests indicated that the computed D P is sensitive to seismic attenuation and can be applied to identify hydrocarbon enrichment in shale gas reservoirs. Obtained D ε reflects the frequency-dependent anisotropy of shale and thus acts as an indicator of bedding fractures. The method was applied to the Longmaxi shale gas reservoir. Results showed that the predicted D P values agreed with the reservoir properties measured in gas-producing boreholes. A higher D P anomaly corresponds to higher TOC content, porosity, and gas content in the target shale gas reservoir. The D ε values estimate the development intensity of bedding fractures, with the results validated by FMI logs and rock physical analyses. Therefore, the proposed dispersion attributes can be used for hydrocarbon identification and bedding fracture detection in organic shale and are thus valuable for shale gas characterization. However, owing to the lack of sophisticated rock physical modeling methods, the influence of organic matter and fractures on the frequency-dependent anisotropy has not been Fig. 17 Map of predicted favorable zones (in red) in the Longmaxi shale gas reservoir considering hydrocarbon enrichment and bedding fractures evaluated by D P and D ε in Fig. 12 fully understood, especially for deep-buried shale. The advancement of rock physics models and laboratory measurements could provide better insight into shale properties and help develop other practical seismic attributes for shale gas characterization. Additionally, vertical fractures were not considered in this study and can be explored in the future using wide-azimuth seismic data. Finally, the structure-related shale properties revealed by the predicted dispersion attributes deserve further investigation based on an improved geological understanding of the local area.