Application of Chemical Sequence Stratigraphy to the Prediction of Shale Gas Sweet Spots in the Wufeng and Lower Longmaxi Formations within the Upper Yangtze Region

: Effective shale gas exploration is hindered by the need for obtaining high-resolution correlations between shale strata and the need for classifying shale facies. To address these issues, chemostratigraphy, sequence stratigraphy, and shale gas geology methods were integrated to develop a new method known as “chemical sequence stratigraphy,” which was successfully applied to the Wufeng–Lower Longmaxi Formations in the upper Yangtze region. Well Huadi 1 was used as a case study, and detailed data were acquired. Multivariate statistical analyses were applied to three deﬁned indices having different genetic signiﬁcance, namely: terrigenous input intensity (TII), authigenic precipitation intensity (API), and organic matter adsorption and reduction intensity (OARI). By analyzing the trends of these three indices, the Wufeng–Lower Longmaxi Formations were divided into ﬁve fourth-order chemical sequences (from bottom to top): LCW, MCL1-1, MCL1-2, MCL1-3, and MCL1-4. The geochemical facies were named and classiﬁed using the chemical sequence stratigraphic framework. The enrichment factor (EF) transformation of elements was conducted to determine whether an element is rich or deﬁcient. The results showed that the favorable geochemical facies in the well were EF-Al deﬁcient, EF-Ca rich, and EF-V rich. The organic matter content and rock brittle strength were then used as chemical parameters, and it was predicted that the LCW and MCL1-1 chemical sequences most likely comprised shale gas sweet spots. This conclusion is consistent with the drilling results and indicates that our proposed method is effective and reliable. This method is further applied to the Changning Shuanghe section, the Shizhu Liutang section, and sections in the Xindi 1 well in the upper Yangtze region. The comparative study of these four sections showed that LCW and MCL1-1 are the key chemical sequences for shale gas exploration and development in the Wufeng–Lower Longmaxi Formations within the Upper Yangtze region.


Introduction
The world is facing huge challenges associated with global warming caused by excessive carbon emissions. Shale gas is one of the cleanest fossil energy sources, and its use is thus significant for reducing such emissions. During the 13th Five-Year Plan (2016-2020) period, China discovered more than 500 billion m 3 of proven geological shale gas field reserves in Changning, Weiyuan, Zhaotong, Fuling, and Weirong [1]. The Wufeng-Lower Longmaxi Formations are composed of deepwater shelf facies, and they are the main gas producing layer within these super-large shale gas fields [1][2][3][4][5][6][7][8]. In recent years, shale gas production has been increasing rapidly, from 108 × 10 8 m 3 in 2018, to 154 × 10 8 m 3 in 2019, and over 200 × 10 8 m 3 in 2020 [8].

Theories and Methods
Sedimentary strata provide a record about the sediment source, sedimentary environment, diagenesis, and chemical element signals. The basic premise of our new method, chemical sequence stratigraphy, is to divide sequences and system tracts using rules It is well known that sea level rise and fall affects changes in the sedimentary environment. In addition, different rock types are formed in different sedimentary environments. These rocks contain various minerals, and the minerals are composed of constituent elements (Table 1); therefore, the elements have genetic significance. For example, Al, K, Ti, and Zr are the most common representative elements of a terrigenous clastic input [26][27][28][29][30], Ca, Mg, Mn, and Ba are the most common representative elements of authigenic precipitation [10,31], and V, U, and Mo are the most common representative elements of reduction environment [32][33][34]. The different chemical elements within strata have particular responses to the relative rise and fall of sea level, and the element combinations therein can be considered as an effective response index. In this study, sequences were divided based on such responses. The geochemical facies were divided based on the chemical sequence stratigraphic framework, and the shale gas sweet spots were ultimately predicted. Table 1. Genetic relationship between elements and minerals (modified from [10,33] The key steps included the following: (1) compiling an optimum geochemical element combination index; (2) establishing a chemical sequence stratigraphy theory model, and dividing the sequence and system tract; (3) classifying and naming geochemical facies and selecting geochemical facies that are favorable for shale gas exploration; (4) evaluating the key chemical parameters of the shale gas sweet spot.

Optimum Geochemical Element Combination Index
The first step used to develop our method was to establish an optimum geochemical element index. According to chemical stratigraphy research, geochemical elements in marine shale are mainly derived from either authigenic mineral precipitation or crystallization, or from terrigenous detrital, biological enrichment, or submarine hydrothermal inputs [9,10,30,[32][33][34]. Of these, submarine hydrothermal inputs are unrelated to changes in sea level, and we thus used three indices to compile the optimum geochemical element combination index: terrigenous input intensity (TII), authigenic precipitation intensity (API), and organic matter adsorption and reduction intensity (OARI). These three indices are described as follows.
(1) Terrigenous input intensity (TII) is a chemical element index that reflects changes in the total terrigenous debris input, and the terrigenous debris input intensity is closely related to the relative rise and fall of regional sea level. With a relative rise in sea level, the area covered by seawater gradually expands, the area of exposed provenance gradually shrinks, the corresponding climate becomes gradually more humid, and the input intensity of terrigenous debris gradually weakens, whereas the opposite occurs with a relative decrease in sea level [10,32]. Because relative sea level rise and fall is regional (and even global), TII was the main index used to analyze chemical sequence stratigraphy.
(2) Authigenic precipitation intensity (API) is a chemical element index that reflects changes in the total amounts of authigenic precipitated or crystalline minerals in the original sedimentary environment. Argillaceous sediments are mainly composed of terrigenous clasts and authigenic precipitated minerals, and their proportional contents have opposite variations; with a relative rise of sea level, the input intensity of terrigenous detritus gradually decreases, while the intensity of authigenic precipitation increases, and the opposite occurs when the relative sea level decreases [10,[32][33][34]. Authigenic precipitates or crystalline minerals are formed in both original sedimentary environments and burial diagenetic environments, and the latter does not reflect relative changes in sea level. Therefore, API was selected as an auxiliary index to analyze chemical sequence stratigraphy.
(3) Organic matter adsorption and reduction intensity (OARI) is a chemical element index that reflects changes in the redox intensity and the total amount of adsorbed elements in organic matter and organic-rich clay. Organic matter is usually enriched and preserved in an anoxic reductive environment, but is easily decomposed by oxidizing, or aerobic bacteria in an oxidizing environment [35]. With the relative rise in sea level, there is an increase in the hypoxic area of deep water and a subsequent gradual increase in organic matter, and the opposite occurs with a relative decrease in relative sea level [10,32,33]. As the elements within organic matter and organic-rich clay can also be derived from terrigenous heavy minerals and carbonate minerals, which cannot completely reflect changes in oxidation reducibility, the OARI was also used as an auxiliary index for chemical sequence stratigraphy in this study.
We also used multivariate statistical analyses that included a correlation analysis, principal component analysis, and cluster analysis to select element combinations from the many major and trace elements in shale formations that had similar genetic significance. These three analysis types were used to determine whether the combinations of elements were highly correlated, had a high load value, and had a small Euclidean distance, respectively [9,32,33,35].

Establishing the Chemical Sequence Stratigraphy Theory Model
The second step in developing our method was to establish a theoretical chemical sequence stratigraphy model that could be used to divide the sequence and system tracts according to the optimized geochemical element combination index. Chemical sequence stratigraphy is consistent with the transgressive-regressive (T-R) model of classical sequence stratigraphy. The basic theory of the T-R model is that taking the lowest point of relative sea level decline as the sequence boundary, the highest point of relative sea level rise as the maximum sea flooding surface, and a sedimentary sequence with a relative sea level rise and fall cycle as a sequence, the sequences can be divided into transgressive system tracts and regressive system tracts [36][37][38]. Chemical sequence stratigraphy also responds to the relative rise and fall of sea level, which is similar to the T-R model, and the sequences division method of chemical sequence stratigraphy is not determined by lithology, but by geochemical element combination indices. In transgressive systems with a rising sea level, the TII index shows a decreasing trend, while the API and OARI indices show increasing trends. In regressive systems with a declining sea level, the TII index shows an increasing trend, while the API and OARI indices show decreasing trends. The sequence boundaries are located at the maximum value of the TII index and the minimum values of the API and OARI indices, and the maximum sea flooding surface is located at the minimum value of the TII index and the maximum values of the API and OARI indices ( Figure 1). rise as the maximum sea flooding surface, and a sedimentary sequence with a relative sea level rise and fall cycle as a sequence, the sequences can be divided into transgressive system tracts and regressive system tracts [36][37][38]. Chemical sequence stratigraphy also responds to the relative rise and fall of sea level, which is similar to the T-R model, and the sequences division method of chemical sequence stratigraphy is not determined by lithology, but by geochemical element combination indices. In transgressive systems with a rising sea level, the TII index shows a decreasing trend, while the API and OARI indices show increasing trends. In regressive systems with a declining sea level, the TII index shows an increasing trend, while the API and OARI indices show decreasing trends. The sequence boundaries are located at the maximum value of the TII index and the minimum values of the API and OARI indices, and the maximum sea flooding surface is located at the minimum value of the TII index and the maximum values of the API and OARI indices ( Figure 1).

Figure 1.
Theoretical chemical sequence stratigraphy model. In transgressive systems, the TII decreases, while the API and OARI increase; in regressive systems, the TII increases, while the API and OARI decrease.

Dividing Geochemical Facies
Large sets of thick shale are generally formed in a deepwater environment, and sedimentary facies are generally stable. We cannot distinguish shale gas favorable facies Theoretical chemical sequence stratigraphy model. In transgressive systems, the TII decreases, while the API and OARI increase; in regressive systems, the TII increases, while the API and OARI decrease.

Dividing Geochemical Facies
Large sets of thick shale are generally formed in a deepwater environment, and sedimentary facies are generally stable. We cannot distinguish shale gas favorable facies from monotonous facies by traditional sedimentary study. Although some sedimentary microfacies can be divided, the results are still far from the requirements for identifying shale gas sweet spots. The third step used to develop our method was to divide and name the geochemical facies. The geochemical facies were classified based on the chemical sequence stratigraphic framework, and these were more precise than sedimentary facies.
Geochemical facies contain main geochemical element combinations that reflect changes in sedimentary environmental conditions [39]. To analyze geochemical facies using the method developed, Al standardization and enrichment factor (EF) transformation of elements are conducted, which are explained as follows. (1) The main aim of the Al standardization of elements is to attenuate the effects of the terrigenous detrital input and amplify the effects of the non-terrigenous input element [40,41], calculated as Al standardized value of the element (X): Al-X = X/Al.
(2) The main aim of the enrichment factor transformation of elements is to clearly determine whether an element is rich or deficient within a sample [33], and the calculation formula is the ratio of the standardized Al value of an element in the shale samples to the standardized Al value of the same element of the world average shale [42][43][44], EF transformative value of the element (X): EF-X = (X/Al) sample /(X/Al) AS , where AS represents "global average shale." When EF-X is equal to 1, the element X of this sample has the same enrichment degree as that of globally average shale; when EF-X is less than 1, X of this sample is less deficient than that of globally average shale; and when EF-X is greater than 1, X of this sample is richer than that of globally average shale [10]. In our method, we divide the value of EF-X into six levels: EF-X < 0.1 is extremely deficient, 0.1 ≤ EF-X < 0.5 is deficient, 0.5 ≤ EF-X < 1 is slightly deficient, 1 ≤ EF-X < 2 is slightly rich, 2 ≤ EF-X < 5 is rich, and EF-X ≥ 5 is extremely rich. In general, a depositional environment with a low terrigenous detrital input intensity, low authigenic precipitation intensity, and a high adsorption and reduction organic matter intensity is favorable for the enrichment of organic matter in shale, and vice versa. In other words, when the TII index is more deficient, the API index is more deficient, and the OARI index is richer, we know that the organic matter is more enriched. Conversely, when the TII index is richer, the API index is richer, and the OARI index is more deficient, the organic matter is more deficient.
To reduce the variables, the most sensitive elements in TII, API, and OARI are further screened out. In general, the most iconic element of TII is Al, the most iconic element of API is Ca, and the most iconic element of OARI is V, U, or Th [9,10,30,[32][33][34]. Organic matter is most enriched in geochemical facies when EF-Al is extremely deficient, EF-Ca is extremely deficient, and EF-V, EF-U, or EF-TH are extremely rich; such facies are the most favorable geochemical facies for shale gas exploration.

Evaluation Key Chemical Parameters
The final step in the model process is to evaluate the key chemical parameters of the shale gas sweet spot based on chemical facies division.
The shale gas sweet spot is a shale segment that has a high gas content, with high productivity [9,45]. The shale gas content has a very close relationship with the total organic carbon (TOC) content, and the brittleness of shale is an important factor affecting shale gas productivity [46,47]. In our model, the TOC is described using the TOC geochemical index, and shale brittleness is described using (SiO 2 + CaO) as a substitute geochemical index. The shale gas potential can be evaluated and the sweet spot predicted based on these two parameters, which are described as follows.
(1) Total organic carbon (TOC) index: Results show that the TOC is positively correlated with shale gas content. In the lower Silurian shale in southern China, when the TOC is >1.0%, the shale gas content is >1.2 m 3 /t; when the TOC is >2.0%, the shale gas content is > 2.0 m 3 /t; and when the TOC is >3.0%, the shale gas content is >4.0 m 3 /t [1]. The shale gas content can be determined by calibrating the TOC measured data, and the its potential can be predicted.
In this study, The shale gas potential is divided into three levels, according to the TOC: when TOC is < 1%, the shale gas potential is average; when 1% ≤ TOC < 3%, the shale gas potential is good; and when the TOC is ≥3%, the shale gas potential is excellent.
(2) (SiO 2 + CaO) index: The matrix permeability of shale is low (generally less than 0.1 × 10 −3 µm 2 ), and the existence of fractures can significantly improve shale gas produc-tivity [48]. Shale segments that have a high brittle mineral content and strong brittleness are conducive to the formation of natural fractures and artificial fracturing [49]. The brittle mineral content of Wufeng-Longmaxi shale in the Changning gas field, the Weiyuan gas field, the Taiyang gas Field, and the Fuling gas field is 42-95%, 42-96%, 55-70%, and 50-80%, respectively [8]. Its average value is therefore approximately 66%, which shows that the brittle mineral content is relatively high. The brittle minerals in shale are non-clay minerals and are mainly silicate and carbonate minerals. The SiO 2 content has a good positive correlation with the silicate mineral content, and the CaO content has a good positive correlation with the carbonate mineral content [8,50]. Therefore, the relative content of (SiO 2 + CaO) can represent the relative contents of silicate and carbonate minerals, and as such, can be used as a substitute index for the brittle mineral content in shale.
In this study, the brittleness of shale is divided into two grades according to the relative content of (SiO 2 + CaO): when (SiO 2 + CaO) < 70%, the shale has low brittleness, and when (SiO 2 + CaO) ≥ 70%, the shale is highly brittle. In conclusion, shale layers that have an excellent potential (TOC ≥ 3%) or a good potential (1% ≤ TOC < 3%) and high brittleness ((SiO 2 + CaO) ≥ 70%) can be predicted, and shale sweet spots can be determined.

Section and Samples
The principles and methods used in chemical sequence stratigraphy were applied to four sections with large spans in the upper Yangtze region, namely the Changning Shuanghe (CS) section, the Shizhu Liutang (SL) section, and sections in the Xindi 1 (XD1) well and the Huadi 1 (HD1) well ( Figure 2). The CS section and SL section are located approximately 348 km apart in the northwest margin of the Dianqianchuane fold belt. HD1 and XD1 are located approximately 340 km apart in the Sichuan basin. In the northeast, the SL section is located approximately 135 km from HD1, and in the southwest, the CS section is located approximately 80 km from XD1. A comparative analysis of the chemical sequence stratigraphy, geochemical facies, and key indices was conducted between the sections, and shale gas sweet spots were predicted in the Wufeng-Longmaxi Formations within the upper Yangtze region. Taking HD1 as an example, the application of chemical sequence stratigraphy in shale gas sweet spot prediction is briefly discussed in the following paragraphs.
HD1 is located in the high part of the Sihaishan-Honghuadian segment within the Huayingshan anticline, which is located in the high and steep fold belt of the eastern Sichuan Basin, where the Triassic is exposed (Figure 2). From bottom to top, the following are revealed: the middle Ordovician Shizipu Formation, the upper Ordovician Baota Formation, the Linxiang Formation, the Wufeng Formation, the Silurian Llanovery Longmaxi Formation, the Shiniulan Formation, the Silurian Wenlock Hanjiadian Formation, the upper Carboniferous Huanglong Formation, the lower Permian Liangshan Formation, the middle Permian Qixia Formation, the Maokou Formation, the upper Permian Longtan Formation, and the Changxing Formation. The Ordovician Wufeng Formation and the Silurian Longmaxi Formation are the main drilling targets.
The Wufeng Formation (O 3 w) has a thickness of 13.74 m, is divided into upper and lower sections, and is in integrated contact with the underlying Linxiang Formation. The Longmaxi Formation (S 1 l) has a thickness of 486.4 m, is divided into upper and lower sections, and is in integrated contact with the underlying the Wufeng Formation and overlying the Shiniulan Formation. Lithologically, the Wufeng Formation mainly comprises calcium-siliceous mudstone, siliceous mudstone, carbonaceous mudstone, and silty mudstone. The lower section of the Longmaxi Formation comprises carbonaceous shale, shale, calcareous mudstone, and interbedded silty mudstone, with unequal thickness that is dominated by black shale. The upper section of the Longmaxi Formation comprises silty carbonaceous mudstone, carbonaceous silty mudstone, silty mudstone, mudstone, and carbonaceous mudstone. With respect to their sedimentary environments, the lower section of the Wufeng Formation is a silica-rich argillaceous deepwater shelf, the upper section of the Wufeng Formation and the lower section of Longmaxi Formation are argillaceous deepwater shelf microfacies, and the upper section of the Longmaxi Formation is a silt-rich argillaceous shallow water shelf microfacies ( Figure 3). HD1 is located in the high part of the Sihaishan-Honghuadian segment within the Huayingshan anticline, which is located in the high and steep fold belt of the eastern Sichuan Basin, where the Triassic is exposed (Figure 2). From bottom to top, the following are revealed: the middle Ordovician Shizipu Formation, the upper Ordovician Baota Formation, the Linxiang Formation, the Wufeng Formation, the Silurian Llanovery Longmaxi Formation, the Shiniulan Formation, the Silurian Wenlock Hanjiadian Formation, the upper Carboniferous Huanglong Formation, the lower Permian Liangshan Formation, the middle Permian Qixia Formation, the Maokou Formation, the upper Permian Longtan Formation, and the Changxing Formation. The Ordovician Wufeng Formation and the Silurian Longmaxi Formation are the main drilling targets.
The Wufeng Formation (O3w) has a thickness of 13.74 m, is divided into upper and lower sections, and is in integrated contact with the underlying Linxiang Formation. The Longmaxi Formation (S1l) has a thickness of 486.4 m, is divided into upper and lower sections, and is in integrated contact with the underlying the Wufeng Formation and overlying the Shiniulan Formation. Lithologically, the Wufeng Formation mainly comprises calcium-siliceous mudstone, siliceous mudstone, carbonaceous mudstone, and silty mudstone. The lower section of the Longmaxi Formation comprises carbonaceous shale, shale, calcareous mudstone, and interbedded silty mudstone, with unequal thickness that is dominated by black shale. The upper section of the Longmaxi Formation comprises silty carbonaceous mudstone, carbonaceous silty mudstone, silty mudstone, mudstone, and carbonaceous mudstone. With respect to their sedimentary environments, the lower section of the Wufeng Formation is a silica-rich argillaceous deepwater  In HD1, the Wufeng-Lower Longmaxi Formations were cored continuously, and 51 samples were collected. The samples from the Wufeng Formation were numbered D1-D14, and those from the Lower Longmaxi Formation were numbered D15-D51 ( Figure  3). In HD1, the Wufeng-Lower Longmaxi Formations were cored continuously, and 51 samples were collected. The samples from the Wufeng Formation were numbered D1-D14, and those from the Lower Longmaxi Formation were numbered D15-D51 (Figure 3).

Data Processing and Analysis
Mineralogical and geochemical tests were conducted on 51 samples obtained from HD1. Mineralogical tests of the samples were conducted using thin section authentication and X-ray diffraction (XRD), and the major and trace elements were determined using a wavelength dispersive X-ray fluorescence (XRF) spectrometer. The analysis was finished in Experimental Center of Research Institute of China University of Geosciences (Beijing).

Data Processing
Data processing included zero value processing and outlier value processing, which are described as follows.
(1) Zero value processing: The zero value in the data is not the real zero value; rather, it is a missing value caused by the content of certain instrument [51]. chemical elements in the sample being less than the minimum detection value of the detection The existence of a zero value has an important influence on evaluating parameters and analyzing the structure of data, and it can destroy the ability to determine a real relationship between components or variables. In particular, it is not possible to conduct a subsequent logarithmic transformation of data [52,53]. Methods used to manage the effects of zero values (missing values) in data mainly include the correlation variable correction method, the maximum expectation algorithm, and the empirical value method. For the data of HD1, the correlation variable correction method was adopted. When the zero value appears in part of the data of element X, another element of the same sample needs to be chosen as an aim element. The selection criteria for the aim element is that it has no zero value and has a high correlation coefficient (R 2 ) with element X. To replace the missing value of element X, an appropriate value that was lower than the lowest detected value was obtained by multiplying the data of the aim element by the given coefficient.
(2) Outlier value processing: Outlier values can significantly affect the results of statistical calculations, especially those based on covariance [54,55]. Multivariate statistical methods are generally used to detect and manage outlier values. It is generally assumed that the values of variables have a normal or lognormal distribution [56]. In this study, no outlier values were found after zero value processing.

Data Analysis and Index Optimization
The following minerals were detected in HD1: clay minerals, quartz, feldspar, calcite, pyrite, and mica. The main elemental oxides detected were Fe 2 O 3 , Al 2 O 3 , K 2 O, Na 2 O, MgO, TiO 2 , ZrO 2 , CaO, SiO 2 , and MnO. The trace elements detected were Zr, V, Cr, Ni, Zn, Cu, Ba, P, As, Pb, Sr, S, and Cl. The correlation analyses of these elements and minerals were conducted. Although SiO 2 had a very good positive correlation with quartz (R 2 = 0.892), it could not be selected as an indicator because of its unclear chemical stratigraphy significance. It was considered possible that the high correlation could be a response to both terrigenous clastic quartz and to chemical and biological quartz, but the chemical stratigraphic significance of the two types of quartz would be quite different. The main terrigenous minerals were feldspar and clay minerals, and these were positively correlated with Al 2 O 3 , K 2 O, TiO 2 , and Cr; for example, the R 2 value between clay minerals and Al 2 O 3 was as high as 0.931. The main authigenic precipitated minerals were calcite, dolomite, and pyrite, and there were obvious positive correlations with CaO, MnO, and P 2 O 5 . For example, the R 2 value of calcite and CaO was as high as 0.897. The R 2 value between TOC and V, Ni, and Zn was above 0.5, and the R 2 value between TOC and V was as high as 0.678. These results show that Al was the typical terrigenous input indicator, Ca was a typical authigenic precipitation indicator, and V was a typical organic matter adsorption indicator. These calculation results are basically consistent with those summarized in previous studies (Table 1). However, to optimize the combination of elements that had the same genetic significance, SPSS statistical software was used to analyze the data, and the following analyses were conducted.
(1) Correlation analysis: As shown in Table 2, the R 2 value of Al 2 O 3 and K 2 O was 0.92, the R 2 value of Al 2 O 3 and TiO 2 was 0.61, and the R 2 value of Cr and TiO 2 was 0.60. Furthermore, the R 2 value of CaO with MgO, MnO, and Sr was 0.53, 0.48, and 0.82, respectively, and the R 2 value of MgO and MnO was as high as 0.89. The R2 value of V with Ba and Cr was 0.91 and 0.81, respectively, and the R 2 value of Ni and Zn reached over 0.5 (0.59), which was very high compared to Ni and other elements. (2) Clustering analysis: Figure 4 shows that the square Euclidean distance (SED) between Al 2 O 3 and K 2 O was approximately 3, between TiO 2 and Cr was approximately 5, and between the combination of Al 2 O 3 , K 2 O, TiO 2 , and Cr was less than 7. The SED of MnO and MgO was 2, the SED between MgO, MnO, and CaO was 4, and the SED between CaO, MgO, MnO, and Sr was 17. The SED of V and Ba was the lowest at 2. The SED of Ni and P was 13, and the overall SED between V, Ba, and Zn was within 17. Generally, element combinations from the many major and trace elements in shale formations that had similar genetic significance had a small Euclidean distance.
(3) Clustering analysis: As evident from Tables 3 and 4 As mentioned above, the same element combinations have the following characteristics: a high R 2 , a small SED value, and a high load value. In summary, the optimal combinations in HD1 associated with the three indices are as follows: (1) for the terrigenous input intensity (TII) index: Al, K, Ti, Cr; (2) for the authigenic precipitation intensity (API) index: Ca, Mg, Mn, Sr; and (3) for the organic matter adsorption and reduction intensity (OARI) index: V, Ni, Ba, Zn.
between Al2O3 and K2O was approximately 3, between TiO2 and Cr was approximately 5, and between the combination of Al2O3, K2O, TiO2, and Cr was less than 7. The SED of MnO and MgO was 2, the SED between MgO, MnO, and CaO was 4, and the SED between CaO, MgO, MnO, and Sr was 17. The SED of V and Ba was the lowest at 2. The SED of Ni and P was 13, and the overall SED between V, Ba, and Zn was within 17. Generally, element combinations from the many major and trace elements in shale formations that had similar genetic significance had a small Euclidean distance.

Chemical Sequence Stratigraphy Division of HD1
Different units of measurement were used in the calculations due to the large differences in the order of magnitude between the major and trace element contents. In this respect, the TII index adopts the original data of the elements, the API adopts the Al standardized value (X/Al) of the elements, and the OARI adopts the original data and an enrichment factor transformative value (EF-X) of the elements. For HD1, the TII index values were the (Al + K + Ti) total content and Cr content, the API index values were the (Ca + Mg + Mn)/Al and Sr/Al contents, and the OARI index values were the (V + Ni + Ba + Zn) total content and the EF-(V + Ni + Ba + Zn).
Based on a comprehensive analysis of the trend in variation of the TII, API, and OARI indices, one chemical sequence (named LCW) was identified within the Wufeng Formation of HD1. Four chemical sequences (named MCL1-1, MCL1-2, MCL1-3, and MCL1-4, from the bottom up) were identified in the Lower Longmaxi Formation of HD1; the respective interfaces of each sequence were named SBW, SBL1, SBL1-2, SBL1-3, and SBL1-4, and the maximum flooding surfaces of each sequence were named mfsW, mfsL1-1, mfsL1-2, mfsL1-3, and mfsL1-4, respectively. The time span of these sequences is less than 1Ma. According to research results on the sequence order [57,58], they are equivalent to fourthorder sequences, and are more highly accurate than the third-order sequences previously identified. The data of each sequence and the key interfaces are shown in Table 5.
According to the histogram for the chemical sequence stratigraphy analysis conducted in HD1 (Figure 5), the terrigenous input intensity (TII) index showed an obvious increasing trend from bottom to top, which reflected a gradual increase in terrigenous input. The TII index increased relatively near each sequence interface, decreased relatively near each maximum sea flooding surface, and underwent a cyclical change from decreasing to increasing in each sequence. This cyclical change represents a response to regional sea level changes, and it can thus be used as the basis for a regional stratigraphic correlation. As previously mentioned, TII was the main index used in chemical sequence stratigraphic division and correlation.  The authigenic precipitation intensity (API) index showed a decreasing trend from bottom to top. Generally, the API index decreased relatively near the interfaces of SBW, SBL1, and SBL1-2, increased relatively near the interfaces of mfsW, mfsL1-1, and mfsL1-2, and showed cyclic changes from increasing to decreasing in the LCW and MCL1-1 chemical sequences. However, the opposite API index changes were seen in the upper MCL1-2, MCL1-3, and MCL1-4, which could be related to the proximity of the The authigenic precipitation intensity (API) index showed a decreasing trend from bottom to top. Generally, the API index decreased relatively near the interfaces of SBW, SBL1, and SBL1-2, increased relatively near the interfaces of mfsW, mfsL1-1, and mfsL1-2, and showed cyclic changes from increasing to decreasing in the LCW and MCL1-1 chemical sequences. However, the opposite API index changes were seen in the upper MCL1-2, MCL1-3, and MCL1-4, which could be related to the proximity of the HD1 to the central Sichuan uplift. When the water is shallow, the API index is affected by injections of fresh water and evaporation. When the sea level falls, seawater salinity increases, and autogenic precipitation is enhanced; when the sea level rises, seawater salinity decreases, and autogenic precipitation is weakened. As previously mentioned, the API was used as an auxiliary index in this study.
The organic matter adsorption and reduction intensity (OARI) index showed a decreasing trend from bottom to top. It decreased relatively near the interfaces of SBW and SBL1, and increased relatively near the mfsW interface. In the LCW chemical sequence, the OARI changed cyclically from increasing to decreasing. However, in the upper MCL1-1, lower MCL1-2, upper MCL1-3, and MCL1-4 chemical sequences, the changes were opposite. This trend was consistent with that of the TII index, and was probably related to shallow water and the influence of terrigenous organic matter input. As mentioned above, the OARI was also employed as an auxiliary index.

Geochemical Facies Division of HD1
For the Wufeng-Lower Longmaxi Formations in HD1, the most iconic elements in the TII, API, and OARI indices were Al, Ca, and V, respectively. Therefore, under the chemical sequence stratigraphic framework, the geochemical facies of the well were named according to the enrichment and deficiency degrees of EF-Al, EF-Ca, and EF-V. According to the classifications used in this study, the geochemical facies of HD1 were determined as follows (Table 6, Figure 6). Table 6. Geochemical facies division of Wufeng-Lower Longmaxi Formation in Well Huadi 1 (SBsequence boundary; msf-maximum sea flooding surface; 0.1 ≤ EF-X < 0.5 is deficient, 0.5 ≤ EF-X < 1 is slightly deficient, 1 ≤ EF-X < 2 is slightly rich, 2 ≤ EF-X < 5 is rich, and the medium means EF value is approximately 1).

Geochemical Facies Number Position
Geochemical Facies Description In the lower part of LCW slightly deficient slightly rich slightly rich 2 In the vicinity of mfsW deficient slightly rich rich 3 From SBL1 to msfL1-1 slightly deficient rich slightly rich 4 In the vicinity of SBL1-2 slightly deficient slightly rich slightly rich 5 From mfsL1-2 to SBL1-3 slightly deficient rich slightly rich 6 From mfsL1-3 to SBL1-4 medium slightly deficient rich 7 In the vicinity of mfsL1-3 slightly rich slightly deficient medium The overall change rule of HD1 is as follows: EF-Al is deficient in the lower and middle parts, rich in the upper part, and is gradually enriched from bottom to top. In contrast, EF-Ca and EF-V gradually decrease from bottom to top.
For HD1, 29 TOC samples from the Wufeng-Lower Longmaxi Formations were tested, all of which were found to have a TOC content greater than 2.0%. The following information was also determined: there was a certain negative correlation between EF-Al and TOC (R 2 = −0.189962), a slightly positive correlation between EF-V and TOC (R 2 = 0.090624), and a minimal correlation between EF-Ca and TOC (R 2 = −0.00144379). According to this relationship, the favorable geochemical facies of shale gas were considered to be EF-Al deficient, EF-Ca rich, and EF-V rich. It can be seen that these favorable geochemical facies are mainly developed in the LCW, MCL1-1, and MCL1-2 chemical sequences (Figure 6), and their measured TOC is higher, with average values of 3.27%, 3.59%, and 3.33%, respectively. The overall change rule of HD1 is as follows: EF-Al is deficient in the lower and middle parts, rich in the upper part, and is gradually enriched from bottom to top. In contrast, EF-Ca and EF-V gradually decrease from bottom to top.
For HD1, 29 TOC samples from the Wufeng-Lower Longmaxi Formations were tested, all of which were found to have a TOC content greater than 2.0%. The following information was also determined: there was a certain negative correlation between EF-Al and TOC (R 2 = −0.189962), a slightly positive correlation between EF-V and TOC (R 2 = 0.090624), and a minimal correlation between EF-Ca and TOC (R 2 = −0.00144379). According to this relationship, the favorable geochemical facies of shale gas were considered to be EF-Al deficient, EF-Ca rich, and EF-V rich. It can be seen that these favorable geochemical facies are mainly developed in the LCW, MCL1-1, and MCL1-2 chemical sequences ( Figure 6), and their measured TOC is higher, with average values of 3.27%, 3.59%, and 3.33%, respectively.

Sweet Spot Prediction based on HD1
The key parameters of TOC and (SiO2 + CaO) were evaluated based on the analysis of favorable geochemical facies. As discussed above, shale sweet spots occur in shale layers with an excellent potential (TOC ≥ 3%) or a good potential (1% ≤ TOC < 3%), and high brittleness ((SiO2 + CaO) ≥ 70%).

Sweet Spot Prediction based on HD1
The key parameters of TOC and (SiO 2 + CaO) were evaluated based on the analysis of favorable geochemical facies. As discussed above, shale sweet spots occur in shale layers with an excellent potential (TOC ≥ 3%) or a good potential (1% ≤ TOC < 3%), and high brittleness ((SiO 2 + CaO) ≥ 70%).
The results of evaluating the shale gas potential in the Wufeng-Lower Longmaxi Formations from HD1 are as follows (Table 7, Figure 7).
Among these good and excellent shale gas potential areas, the lower and middle parts of the LCW chemical sequence, the vicinity of SBL1-2, and the middle and lower parts of the MCL1-2 chemical sequence all contain high brittle shale (Figure 7), and are thus predicted to be the shale gas sweet spots in HD1. Table 7. The shale gas potential of Wufeng-Lower Longmaxi Formation in Well Huadi 1. (SBsequence boundary; msf-maximum sea flooding surface; the shale gas potential is average when TOC < 1%, the potential is good when 1% ≤ TOC < 3%, and the potential is excellent when TOC ≥ 3%).

Number
Position Shale Gas Potential The results of evaluating the shale gas potential in the Wufeng-Lower Longmaxi Formations from HD1 are as follows (Table 7, Figure 7). Table 7. The shale gas potential of Wufeng-Lower Longmaxi Formation in Well Huadi 1. (SB-sequence boundary; msf-maximum sea flooding surface; the shale gas potential is average when TOC < 1%, the potential is good when 1% ≤ TOC < 3%, and the potential is excellent when TOC ≥ 3%).

Number
Position Shale Gas Potential in the vicinity of mfsL1-4 average Figure 7. Histogram of shale gas sweet spot prediction results for the Wufeng-Lower Longmaxi Formations within Well Huadi 1. "The main production layer" and "the biggest total gas content layer" [59] are basically consistent with the sweet spots (yellow) identified in this study).
The wellhead blowout pipeline of HD1 reportedly successfully induced ignition, with a flame height of 7-8 m, and the gas-producing layers are considered to be "the main production layer" and "the biggest total gas content layer" [59]. "The main production layer" is consistent with the shale gas sweet spot identified in the middle and lower part of LCW chemical sequence, and "the biggest total gas content layer" is basically consistent with the sweet spot identified near SBL1-2 and the middle and lower parts of the MCL1-2 chemical sequence (Figure 7). Therefore, these results prove that the research method presented here is feasible.

Comparative Study of the Upper Yangtze
For the fine isochronal division of strata over a large area, the chemical sequence stratigraphy research method was applied to the Changning Shuanghe (CS) section, the Shizhu Liutang (SL) section, and sections in the Xindi 1 (XD1) well and the Huadi 1 (HD1) well in the upper Yangtze region. The depositional background of the Wufeng-Lower Longmaxi Formations in the upper Yangtze region is a craton platform covered by vast amounts of seawater that lacks a terrigenous coarse clastic supply. High-quality dark mudstone is mainly distributed in the Shizhu-Yibin area, where the maximum cumulative thickness reaches over 120 m, but is generally greater than 40 m [60]. The SL section is located in the thickest depositional center of the high-quality dark mudstone, the HD1 well is located near the Huayingshan fault structural belt where the high-quality dark mudstone is the thinnest, and the CS section and the XD1 well are located between them ( Figure 8). The chemical sequence stratigraphy research results for the XD1, CS, and SL sections have been previously published [31,61,62]. Such results were used to conduct comparative research for the HD1 well-XD1 well-CS section-SL section, which runs from a point near the uplift to the deposition center, with the aim of analyzing the change rules for the chemical sequences, the geochemical facies, and the shale gas potentials.  . The SL section is located in the thickest depositional center of the high-quality dark mudstone, the HD1 well is located near the Huayingshan fault structural belt where the high-quality dark mudstone is the thinnest, and the CS section and the XD1 well are located between them.
As seen in the table providing the statistics of chemical sequence thicknesses (Table  8) and the joint comparison between the chemical sequence stratigraphy and geochemical facies (Figure 9), five fourth-order chemical sequences (LCW, MCL1-1, MCL1-2, MCL1-3, and MCL1-4, from the bottom up) were identified in each section over a large area. The same sequences in four sections show obvious differences in the thickness changes, and these are presumed to be mainly related to the terrigenous input and authigenic precipitation and tectonic subsidence and later reformation. For example, the thickness of LCW  Wang et al., 2017). The SL section is located in the thickest depositional center of the high-quality dark mudstone, the HD1 well is located near the Huayingshan fault structural belt where the high-quality dark mudstone is the thinnest, and the CS section and the XD1 well are located between them.
As seen in the table providing the statistics of chemical sequence thicknesses (Table 8) and the joint comparison between the chemical sequence stratigraphy and geochemical facies (Figure 9), five fourth-order chemical sequences (LCW, MCL1-1, MCL1-2, MCL1-3, and MCL1-4, from the bottom up) were identified in each section over a large area. The same sequences in four sections show obvious differences in the thickness changes, and these are presumed to be mainly related to the terrigenous input and authigenic precipitation and tectonic subsidence and later reformation. For example, the thickness of LCW undergoes a gradually increasing trend from the uplift to the deposition center. The thickness of MCL1-1 changes from thickening to thinning, that of MCL1-2 changes from thickening to thinning and then thickening, that of MCL1-3 changes from thickening to thinning, and that of MCL1-4 changes from thickening to thinning. In general, the changes in the thicknesses of each chemical sequence in HD1, which is the closest to the Central Sichuan Uplift, are relatively small.   There are similar characteristics in the geochemical facies within the same chemical sequence of each section. LCW and MCL1-1 are dominated by EF-Al deficient, EF-Ca rich, and EF-V-rich chemical facies, which are shale gas favorable chemical facies. MCL1-2 and MCL1-3 are dominated by EF-Al slightly deficient, EF-Ca slightly rich, and EF-V slightly rich geochemical facies, and MCL1-4 is dominated by EF-Al slightly rich, EF-Ca slightly deficient, and EF-V slightly rich or slightly deficient geochemical facies. Vertically, from bottom to top, EF-Al changes from deficient to rich, and EF-Ca and EF-V changes from rich to deficient. There is obvious EF-Al enrichment and EF-V depletion near each chemical sequence interface, while EF-Al depletion and EF-V enrichment appear near each maximum flood surface (Figure 9).
The joint comparison figure shows certain differences in the shale gas potential and the development of the sweet spot in each section ( Figure 10). The SL section has the thickest and best shale gas potential, followed by the CS section. XD1 is thinner, and HD1 is the thinnest, which shows that the shale gas potential decreases from the deposition center to the uplift. For the shale gas sweet spot, LCW and MCL1-1 are the main development chemical sequences, whereas the numbers of shale gas sweet spots in MCL1-2 and MCL1-3 decrease. The MCL1-4 chemical sequence shows an underdeveloped sweet spot.   In conclusion, LCW and MCL1-1 are the key chemical sequences for shale gas exploration and development in the Wufeng-Lower Longmaxi Formations within the Upper Yangtze region.

Conclusions
The new chemical sequence stratigraphic research method developed in this study is an effective tool for finely dividing thick mud shale strata and making associated comparisons.