Geochemical Characteristics and Chemostratigraphic Analysis of Wufeng and Lower Longmaxi Shales, Southwest China

: The demand for shale gas has propelled researchers to focus on precise and high-resolution stratigraphic divisions for homogeneous shales, of which the late-Ordovician Wufeng (O 3 w) and the early-Silurian Longmaxi (S 1 l) formations in southwest China are two of the best candidates for shale gas exploration in China. However, systematic chemostratigraphic work for these strata is still sparse, and the existing chemostratigraphic work either lack representativeness in terms of the proxies used or are subjective during their division procedures. Thus, automatic division process based on multi proxies and an objective statistical technique was applied to establish a quantitative, high-resolution, and robust chemostratigraphic scheme for the Wufeng and lower Longmaxi shales. The geochemical analysis unveils that the Wufeng and Lower Longmaxi shales show prominent heterogeneities in terrigenous inputs, redox conditions, and paleoproductivity, enabling the potential application of chemostratigraphy to these strata. Based on these heterogeneities, the chemostratigraphic scheme for the Wufeng and Lower Longmaxi shales has been established, and the whole strata could be divided into 13 chemozones using constrained clustering analysis. The chemostratigraphic scheme could not only be comparable to the regional sequence stratigraphic scheme but also more objective and higher-resolution. The high TOC content and brittle minerals within chemozone C1 makes it the most preferable layer for shale gas exploration and development. This research gives a systematic chemostratigraphic analysis on Wufeng and Lower Longmaxi shales, which testiﬁes the feasibility and potential of usage of chemostratigraphy for Chinese shale gas exploration and development.


Introduction
Increasing energy demands have motivated interests in unconventional shale gas, and North America, especially the United States, has gradually initiated a worldwide shale gas revolution. Consequently, characterizing shales has become a constant research focus both in scientific studies (e.g., [1]) and industry production (e.g., [2]), in which the precise and high-resolution division and correlation of mudstone successions has become one of the most critical steps [3,4]. However, the macroscopic homogeneity, mineralogic complexity, and limited availability of biostratigraphic constraints of such rocks make it difficult from the perspective of traditional approaches [2,[5][6][7]. Thus, new methodologies are needed to resolve it.
Chemostratigraphy, which utilizes chemical fingerprints or variations of geochemical traits in sedimentary successions for stratigraphic correlation [8,9], could be such a technique. This is because shales are found to have significant heterogeneities from a geochemistry perspective (e.g., [3,10]). An increasing number of researchers are thus attempting to employ the chemostratigraphic technique to establish the stratigraphic correlation scheme (e.g., [6,7,[9][10][11][12]). Other than that, this new technique allows better discrimination of rock Sichuan Basin (Figure 1). The depositions with coarse grains often lead to biased chemostratigraphic results. Thus, in this study, the Upper Longmaxi is ignored due to its high content of arenaceous compositions. The Wufeng Fm. comprises black siliceous and calcareous shales, and contains abundant Brachiopoda and Radiolaria with a thickness of approximately 10 m. The Lower Longmaxi Fm., with a thickness of over 90 m, is mainly characterized by black arenaceous and carbonaceous shales, while on the top of the Lower Longmaxi, arenaceous components are increased dramatically.

Materials
The Wufeng and the Lower Longmaxi Fm. were well preserved and outcropped on the Liutang section. A total of 192 samples, with an average sampling step of approximately 0.5 m, were collected from the bottom of the Wufeng Fm. to the Lower Longmaxi Fm. Thin sections were made for every sample, and they then were identified to ensure that the samples were fresh and did not suffer from any heavy chemical alterations. Then, all qualified sampled were ground into powder using an agate grinder and subjected to geochemical analysis. The Liutang section, as one of the most typical profiles for studying the Wufeng and Longmaxi Fms [23], is located in Shizhu County, Chongqing Municipality, and is near the Jiaoshiba area which is the largest shale gas field in China. Geologically, Liutang is located in the southern part of the Middle-Upper Yangtze Craton and is adjacent to the southeast Sichuan Basin (Figure 1). The depositions with coarse grains often lead to biased chemostratigraphic results. Thus, in this study, the Upper Longmaxi is ignored due to its high content of arenaceous compositions. The Wufeng Fm. comprises black siliceous and calcareous shales, and contains abundant Brachiopoda and Radiolaria with a thickness of approximately 10 m. The Lower Longmaxi Fm., with a thickness of over 90 m, is mainly characterized by black arenaceous and carbonaceous shales, while on the top of the Lower Longmaxi, arenaceous components are increased dramatically.

Materials
The Wufeng and the Lower Longmaxi Fm. were well preserved and outcropped on the Liutang section. A total of 192 samples, with an average sampling step of approximately 0.5 m, were collected from the bottom of the Wufeng Fm. to the Lower Longmaxi Fm. Thin sections were made for every sample, and they then were identified to ensure that the samples were fresh and did not suffer from any heavy chemical alterations. Then, all qualified sampled were ground into powder using an agate grinder and subjected to geochemical analysis.

Geochemical Analysis Method
The Wavelength Dispersive X-ray Fluorescence Spectrometer (XRF) in the Micro Structure Analytical Laboratory of Beijing was used to acquire the contents of the major and trace elements. The major elements, such as Al, Si, Ca, and Ti, were measured as a percentage of metallic-oxides, and the trace elements, such as V, Cr, Rb, Y were measured in mg/kg (ppm). Calibration was performed with 25 rock reference materials (GBW07101-07125) certified by the Chinese Academy of Geological Sciences. The total organic carbon (TOC) was measured in the laboratory at the Academy of Science, China University of Geosciences (Beijing) using a Total Carbon and Total Nitrogen Analyzer of Jena TOC/TN Multi NC2100S (1300128S), and was scaled with % w/w. Prior to the measurement, the samples were ground into fine powder and then treated with nitric acid (2 N or 2 mol/L) to remove inorganic carbon. The datasets for the geochemical measurements are documented in Supplementary Material Table S1. Average Shale (AS) [24] was used for analytical control and data comparison purposes.
In order to unveil the geochemical characteristics of the succession so as to select representative geochemical proxies used in the chemostratigraphic process, the geochemical analysis was implemented; during this procedure, several geochemical transformations should be calculated. The Al-normalized values (Al-normalized value = element/Al sample ) of elements were firstly calculated to remove the effects of terrigenous inputs ( [25] Calvert and Pederson, 1993). Then, the enrichment factor value EF element was calculated based on Average Shales (AS; [24]) to assess the enrichment level of the elements, and the calculating formula is as follows: EF = (Element/Al) sample ÷ (Element/Al) AS [25]. Finally, the value of excess Si (Si exc ) was calculated using the formula Si exc = [(SiO 2 /Al 2 O 3 ) sample − (SiO 2 /Al 2 O 3 ) AS ] × Al 2 O 3sample (AS = Average Shales). Through this process, the terrigenous inputs were deducted [26].

Scanning Electron Microscopy (SEM)
Scanning electron microscopy (SEM) is an observation method. It uses a narrow focused high-energy electron beam to scan the sample, excites various physical information through the interaction between the beam and the material, and collects, amplifies, and reimages this information to achieve the purpose of characterizing the microscopic morphology of the material [27]. The SEM were carried out on a JEOL JSM-6700F SEM at KYKY Technology Co., LTD, Beijing, China, to investigate the morphology of framboid pyrites.

Multivariate Statistics
After the selection of geochemical proxies, the next most important step for chemostratigraphy is finding a reasonable method for determining the zonation of the geochemical data that results in a reliable and interpretive stratigraphic scheme for the succession of interest [6,[10][11][12][13]18,[28][29][30]. Traditionally, the zonation of geochemical data was conducted in an artificial way that was affected by the experience and knowledge of operators. Consequently, automatic zonation by software programs, especially using the constrained clustering technique, known as hierarchical cluster analysis, was introduced [11,31].
Constrained clustering analysis is a multivariate statistical approach that sorts the hyperdimensional clusters (or samples) based on the similarity coefficient between two neighboring clusters without disordering their original sequence [6]. The similarity refers to the distance between two clusters. In consideration of possible correlations existing among variables of rock samples, the 'Euclidian Distance' method is the preferred way of measuring the distance [6,10]. Thus, the distance between two clusters has been computed based on the length of the straight line drawn between two clusters. The working principle of 'Euclidian Distance' is that it will measure only the direct distance between two centroids. Consequently, the 'K-means method' is thought to be more appropriate. In this method, each cluster is represented by the center of the cluster. In order to maximize the inter-cluster difference and minimize the intra-cluster difference, the 'Minimum Variance' or so-called 'Ward's Method [32] is the optimal one to calculate the centroids of newly-merged clusters in geological statistical applications [10,33,34]. This method makes whichever merger of two clusters, resulting in a minimum increase in total within-cluster variance [33]. The result of constrained clustering analysis is visualized in the form of a dendrogram, in which the geological samples are arranged according to the stratigraphic depth (in meters) [35]. Furthermore, the ultimate number of clusters depends on selecting the cut-off values in the dendrogram [11,36]. The broken-stick model (abbreviated as the bstick model) proposed by Bennett [36] is an efficient way to estimate the statistically significant number of zones in the dendrograms. We used this method to estimate the optimal number of chemozones. The constrained clustering analysis and the bstick model were realized using the Rioja package under CRAN R software.
Prior to multivariate statistical analysis, an EDA (Exploratory Data Analysis, EDA) analysis of raw data was performed in order to identify outliers (<10%) and zeros (<10%), and these were replaced by two-thirds of the determination limits and the corresponding statistical medians, respectively. Notably, the raw datasets, whether trace elements or major elements, were expressed as ppm to conduct the EDA process more conveniently and accurately (e.g., [11]).

Geochemistry
On the Liutang section, Al 2 O 3 , SiO 2 , and CaO were three of the most abundant oxides, their contents are 2.93~20.0%, 32.4~91.9%, and 0.02~17.3%, respectively, with average values of 14.4%, 66.7%, and 0.45%. The contents of Al 2 O 3 and SiO 2 were close to the Average Shale (AS, Al 2 O 3 = 16.6%, and SiO 2 = 58.9%), while the content of CaO was far lower than the Average Shale (AS, CaO = 2.24%). By comparison, it can be deduced that the content of SiO 2 gradually increased from the dark-color shale to light-color shale to arenaceous shale, with the average value of 58.9%, 65.5%, and 64.9%, respectively, while the content of CaO in those rocks decreased from the bottom to the upper profile, with average values of 2.24%, 0.32%, and 0.15%, respectively (see Table 1 for details). Within Wufeng and the Lower Longmaxi Fm., the contents of the terrigenous elements, Rb, Y, and Ti, were generally equivalent to the Average Shale (AS) as reflected by the Alnormalized value and the enrichment factor (EF). Their average EF values were 0.99, 0.92, and 0.96, respectively, while, in the interval of the dark-color shales, their contents showed prominent fluctuating characteristics with EF values of 0.43~2.17, 0.39~3.94, and 0.66~1.78, respectively (see Table 1 for details).
The V/Cr ratio in the Liutang section, reflecting the redox conditions of the bottom water [37], was between 1.05 and 9.28. Moreover, the V/Cr value significantly decreased from the bottom dark-color shales to the middle shales and then to the top argillaceous shales with average values of 4.43, 2.02, and 1.56, respectively. The statistical analysis of the TOC, reflecting total organic carbon within the strata, and the paleo-oceanic productivity indicator, excess Si (Si exc ), show that both the TOC and Si exc decreased from bottom to top. From the dark-color shale to shale to the argillaceous shale, the average TOC and the average Si exc were 3.21%~1.68%~0.31% and 39.44~11.05~5.84, respectively (see Table 1 for details).
The broken-stick model (abbreviated as the bstick model) proposed by Bennett [36] is an efficient way to estimate the statistically significant number of zones in the dendrograms. We used this method to estimate the optimal number of chemozones, which was realized using the bstick.chclust function within the R. As depicted in Figure 2, the numbers 3, 11, and 13 are all statistically significant, and of course, the numbers 2, 4, 6, 8, 10, and 11 are also significant numbers of chemozones. However, too few chemozones (or clusters) may not reveal clear regional patterns [11] and also cannot illustrate the high-resolution features of chemostratigraphy. Thus, we firstly selected the Euclidean distance of 4803 in the dendrogram as the cutoff point, then three chemozones were gained (that is, the three first-level chemozones, C1, C2, and C3), and they were thought to represent the regional stratigraphic pattern, and were regarded as first-level chemozones. Then, to better illustrate the geologically local features of the strata, we moved the cutoff values into 3720, then eleven chemozones (C1-1, C1-2, C1-3, C1-4, C2-1, C2-2, C2-3, C2-4, C2-5, C2-6, C3) were gained, as shown in Figure 4, and designated as the second-level chemozones. Notably, the boundaries between these chemozones correlate well to the changing points of the proxy's profiles ( Figure 4). Nevertheless, some local geochemical turning points ( Figure 4) on profiles of SiO2/(SiO2 + 5Al2O3), V/Cr, TOC, and Siexc were not wholly unveiled, so we shifted the cutoff value into 2573, then four third-level chemozones were obtained (C2-1A, C2-1B, C2-2A, and C2-2B). The boundary between C2-2A and C2-2B was especially prominent on profiles of proxies ( Figure 4). Then, the final number of chemozones was thirteen ( Figure 4), which is statistically significant, just as shown in the bstick model ( Figure 2). On one hand, this hierarchical method unveiled the regional pattern of the strata, and on the other hand, dug into the local features of the strata.     The broken-stick model (abbreviated as the bstick model) proposed by Bennett [36] is an efficient way to estimate the statistically significant number of zones in the dendrograms. We used this method to estimate the optimal number of chemozones, which was realized using the bstick.chclust function within the R. As depicted in Figures 2, the numbers 3, 11, and 13 are all statistically significant, and of course, the numbers 2, 4, 6, 8, 10, and 11 are also significant numbers of chemozones. However, too few chemozones (or clusters) may not reveal clear regional patterns [11] and also cannot illustrate the highresolution features of chemostratigraphy. Thus, we firstly selected the Euclidean distance of 4803 in the dendrogram as the cutoff point, then three chemozones were gained (that is, the three first-level chemozones, C1, C2, and C3), and they were thought to represent the regional stratigraphic pattern, and were regarded as first-level chemozones. Then, to better illustrate the geologically local features of the strata, we moved the cutoff values into 3720, then eleven chemozones (C1-1, C1-2, C1-3, C1-4, C2-1, C2-2, C2-3, C2-4, C2-5, C2-6, C3) were gained, as shown in Figure 4, and designated as the second-level chemozones. Notably, the boundaries between these chemozones correlate well to the changing points of the proxy's profiles ( Figure 4). Nevertheless, some local geochemical turning points ( Figure 4) on profiles of SiO2/(SiO2 + 5Al2O3), V/Cr, TOC, and Siexc were not wholly unveiled, so we shifted the cutoff value into 2573, then four third-level chemozones were obtained (C2-1A, C2-1B, C2-2A, and C2-2B). The boundary between C2-2A and C2-2B was especially prominent on profiles of proxies ( Figure 4). Then, the final number of chemozones was thirteen ( Figure 4), which is statistically significant, just as shown in the bstick model ( Figure 2). On one hand, this hierarchical method unveiled the regional pattern of the strata, and on the other hand, dug into the local features of the strata.

Paleoproductivity
Abundant existing research on the Wufeng and Longmaxi shales have demonstrated that the upper Wufeng and the lower Longmaxi succession are dominated by siliceous shales. Thus, Si could be an optimal candidate for recovering paleoproductivity. Generally, the silica in the strata can be contributed by authigenic (clay-mineral-transformed origin), hydrothermal, terrigenous, and biogenic influxes.
The Si-Al crossplot was used to estimate the clay-mineral-transformed silica. As shown in Figure 5a, almost all samples are plotted above the illite Si/Al line, especially for the samples within Wufeng and Lower Longmax Fms. (83.5-101 m), indicating that silica in the strata should not be attributed to the mineral transformation from smectite into illite. Moreover, the strong negative correlation between Al and Si (−0.66, Figure 5a) also implies a biogenic origin of silica. The Al-Fe-Mn ternary diagram is useful in assessing the difference in chemical composition between hydrothermal and non-hydrothermal siliceous rock [38]. As shown in Figure 5b, almost all samples are plotted at the biogenic area, indicating that silica in the Wufeng and the Lower Longmaxi Fms. has not been affected  Thus, Si could be an optimal candidate for recovering paleoproductivity. Generally, the silica in the strata can be contributed by authigenic (clay-mineral-transformed origin), hydrothermal, terrigenous, and biogenic influxes.
The Si-Al crossplot was used to estimate the clay-mineral-transformed silica. As shown in Figure 5a, almost all samples are plotted above the illite Si/Al line, especially for the samples within Wufeng and Lower Longmax Fms. (83.5-101 m), indicating that silica in the strata should not be attributed to the mineral transformation from smectite into illite. Moreover, the strong negative correlation between Al and Si (−0.66, Figure 5a) also implies a biogenic origin of silica. The Al-Fe-Mn ternary diagram is useful in assessing the difference in chemical composition between hydrothermal and non-hydrothermal siliceous rock [38]. As shown in Figure 5b, almost all samples are plotted at the biogenic area, indicating that silica in the Wufeng and the Lower Longmaxi Fms. has not been affected by hydrothermal influxes [39]. Thus, the diagenetic remobilization of silica has not occurred. Element Zr is typically associated with heavy mineral zircon and can be used as terrigenous proxy. Thus, a Si (%)-Zr (ppm) cross-plot can be used to determine whether the Si in the strata is contributed to by terrigenous or biogenic inputs. A positive correlation between Zr and Si indicates the terrigenous contribution of SiO2, while a negative between them implies biogenic contribution. As shown in Figure 5c, the Si (%)-Zr (ppm) cross-plot shows that the silica within Wufeng Fm. and bottom of the Lower Longmaxi Fm. belong to the biogenic origin, which could be validated by the appearance of siliceous organisms within the strata (see Figure 6j-l). Nevertheless, the silica from the upper part of Lower Longmaxi Fm. is characterized by terrigenous inputs, which is consistent with previous work on geochemical characteristics of Wufeng and Longmaxi Fms. [40,41].
Excess Si (Siexc) is the fraction that deducts the terrigenous inputs [26]. It can be calculated using the formula Siexc = [(SiO2/Al2O3)sample − (SiO2/Al2O3)AS] × Al2O3sample (AS = Average Shales). Moreover, TOC content is the direct reflection of paleoproductivity; high TOC indicates high productivity level. As shown in Figure 5, Siexc and TOC show the same trend ( Figure 5d). Thus, after excluding the possibility of authigenic and hydrothermal contribution, Siexc could be an efficient proxy of paleoproductivity.
Generally, the high Siexc value means a high paleoproductivity level. As shown in Figure 2, the value of Siexc is high at intervals of 83.5~101 m and 62~73 m, indicating that the Wufeng Fm. and Lower Longmaxi Fm. at the corresponding depths have a higher level of paleoproductivity. The siliceous Radiolaria and Spongy spicule (Figure 6j,l) appearing at Wufeng and the bottom of the Lower Longmaxi may be the reason for the increase in Siexc. However, the value of Siexc in the Guanyinqiao layer (GYQ, 91.5~93 m) is significantly lower than in the adjacent layers, indicating that the number of siliceous organisms in this interval has dramatically decreased, which is mainly because the GYQ layer belongs to calcareous submarine sediments [42]. At the intervals of 73~83.5 m and 16~62 m, the value of Siexc is extremely low, indicating a normal-to-low level of paleoproductivity. At the intervals of 0~16 m, the value of Siexc repeatedly fluctuates, indicating that the level of paleoproductivity in the strata has changed frequently at the topmost level of Lower Longmaxi. Element Zr is typically associated with heavy mineral zircon and can be used as terrigenous proxy. Thus, a Si (%)-Zr (ppm) cross-plot can be used to determine whether the Si in the strata is contributed to by terrigenous or biogenic inputs. A positive correlation between Zr and Si indicates the terrigenous contribution of SiO 2 , while a negative between them implies biogenic contribution. As shown in Figure 5c, the Si (%)-Zr (ppm) cross-plot shows that the silica within Wufeng Fm. and bottom of the Lower Longmaxi Fm. belong to the biogenic origin, which could be validated by the appearance of siliceous organisms within the strata (see Figure 6j-l). Nevertheless, the silica from the upper part of Lower Longmaxi Fm. is characterized by terrigenous inputs, which is consistent with previous work on geochemical characteristics of Wufeng and Longmaxi Fms. [40,41].
Excess Si (Si exc ) is the fraction that deducts the terrigenous inputs [26]. It can be calculated using the formula Si exc = [(SiO 2 /Al 2 O 3 )sample − (SiO 2 /Al 2 O 3 ) AS ] × Al 2 O 3sample (AS = Average Shales). Moreover, TOC content is the direct reflection of paleoproductivity; high TOC indicates high productivity level. As shown in Figure 5, S iexc and TOC show the same trend (Figure 5d). Thus, after excluding the possibility of authigenic and hydrothermal contribution, Si exc could be an efficient proxy of paleoproductivity.
Generally, the high Si exc value means a high paleoproductivity level. As shown in Figure 2, the value of Si exc is high at intervals of 83.5~101 m and 62~73 m, indicating that the Wufeng Fm. and Lower Longmaxi Fm. at the corresponding depths have a higher level of paleoproductivity. The siliceous Radiolaria and Spongy spicule (Figure 6j,l) appearing at Wufeng and the bottom of the Lower Longmaxi may be the reason for the increase in Si exc . However, the value of Si exc in the Guanyinqiao layer (GYQ, 91.5~93 m) is significantly lower than in the adjacent layers, indicating that the number of siliceous organisms in this interval has dramatically decreased, which is mainly because the GYQ layer belongs to calcareous submarine sediments [42]. At the intervals of 73~83.5 m and 16~62 m, the value of Si exc is extremely low, indicating a normal-to-low level of paleoproductivity. At the intervals of 0~16 m, the value of Si exc repeatedly fluctuates, indicating that the level of paleoproductivity in the strata has changed frequently at the topmost level of Lower Longmaxi.

Terrigenous Inputs
Al2O3, SiO2, and CaO are three of the most abundant components within shale rocks, their proportion decides the rock types and properties. Al2O3 represents terrigenous inputs [25,43], SiO2 in the strata, as discussed above, is mainly derived from biogenic sources, while CaO is derived from both terrigenous inputs and authigenic oceanic sources [43,44]. In consideration of the low content of Al2O3 compared with SiO2, and to facilitate the comparison with the average shales, five multiples of Al2O3 were used as one member in the 5Al2O3-CaO-SiO2 ternary diagram ( Figure 5e). As shown in Figure 5e, only one sample (sample LT-18D) was plotted on the carbonate dilution line. Thus, the ratio SiO2/(SiO2+5Al2O3) can be further used to ensure whether there are aluminous or siliceous depositions for the remaining samples. Generally, Si/Al is the proxy of biogenic quartz [12,28]. As shown in Figure 5f, a strong negative correlation (−0.97) between Al2O3 (%) and SiO2/(SiO2 + 5Al2O3) emerged, indicating they were mutually diluted between Al2O3 and SiO2. Thus, SiO2/(SiO2 + 5Al2O3) could also be used as a biogenic vs. terrestrial proxy, and the lower ratio indicates more terrigenous inputs and vice versa.
The sample LT-18D represents the lithology of the Guanyinqiao (GYQ) layer or Top Wufeng Fm. As shown in Figure 5e, this sample is diluted by the carbonate, indicating a marine origin of the CaO within the GYQ layer. The numerous Brachiopoda found in the strata reveals that the CaO in the GYQ layer is originated from biogenic sources, while it has terrigenous contributions for other layers. When leaving out sample LT-18D, almost all samples show no effects of carbonate dilution and are similar to Average Shales. In detail, the curve of SiO2/(SiO2+5Al2O3) (Figure 4)  Generally, Rb, Ti, and Y represent terrigenous inputs [43]. Thus, the profiles of Alnormalized values of these elements could also be used to unveil the terrigenous characteristics of the strata. As shown in Figure 7

Terrigenous Inputs
Al 2 O 3 , SiO 2 , and CaO are three of the most abundant components within shale rocks, their proportion decides the rock types and properties. Al 2 O 3 represents terrigenous inputs [25,43], SiO 2 in the strata, as discussed above, is mainly derived from biogenic sources, while CaO is derived from both terrigenous inputs and authigenic oceanic sources [43,44]. In consideration of the low content of Al 2 O 3 compared with SiO 2 , and to facilitate the comparison with the average shales, five multiples of Al 2 O 3 were used as one member in the 5Al 2 O 3 -CaO-SiO 2 ternary diagram ( Figure 5e). As shown in Figure 5e, only one sample (sample LT-18D) was plotted on the carbonate dilution line. Thus, the ratio SiO 2 /(SiO 2 + 5Al 2 O 3 ) can be further used to ensure whether there are aluminous or siliceous depositions for the remaining samples. Generally, Si/Al is the proxy of biogenic quartz [12,28]. As shown in Figure 5f, a strong negative correlation (−0.97) between Al 2 O 3 (%) and SiO 2 /(SiO 2 + 5Al 2 O 3 ) emerged, indicating they were mutually diluted between Al 2 O 3 and SiO 2 . Thus, SiO 2 /(SiO 2 + 5Al 2 O 3 ) could also be used as a biogenic vs. terrestrial proxy, and the lower ratio indicates more terrigenous inputs and vice versa.
The sample LT-18D represents the lithology of the Guanyinqiao (GYQ) layer or Top Wufeng Fm. As shown in Figure 5e, this sample is diluted by the carbonate, indicating a marine origin of the CaO within the GYQ layer. The numerous Brachiopoda found in the strata reveals that the CaO in the GYQ layer is originated from biogenic sources, while it has terrigenous contributions for other layers. When leaving out sample LT-18D, almost all samples show no effects of carbonate dilution and are similar to Average Shales. In detail, the curve of SiO 2 /(SiO 2 + 5Al 2  Generally, Rb, Ti, and Y represent terrigenous inputs [43]. Thus, the profiles of Al-normalized values of these elements could also be used to unveil the terrigenous characteristics of the strata. As shown in Figure 7, the Al-normalized values of Rb, Ti, and Y show similar changing trends. In the interval of 83.5~101 m, the curves of Rb/Al, Ti/Al, and Y/Al show dramatic fluctuation, indicating multiple sources and slow depositional rate in Wufeng and the Bottom Longmaxi Fm. The abundant occurrence of fossils, such as Brachiopods and Radiolaria, within the strata indicate that this fluctuation may be closely related to pulsed biological contribution. In the interval of 83.5~16 m, Rb/Al, Ti/Al, and Y/Al show a slight magnitude of fluctuation indicating that the strata are characterized by a mono provenance and high depositional rate. This characterization is the geochemical response to the fast contemporaneous convergence between 'Dianqiangui' and 'Cathaysia' and a simultaneously fast global regressive event [45,46]. The continuous uplifting of the ancient lands and the consistent warming of climate conditions [47,48] contributed abundant deposition materials to the basin during the deposition of middle Lower Longmaxi Fm. The global regressive event could account for the decline in regional sea levels and increase in the accumulation rate of the strata. Within the interval of 0~16m, Rb/Al, Ti/Al, and Y/Al fluctuate dramatically (Figure 7), indicating that multiple provenance areas contribute deposition materials to the basin simultaneously. Combined with the coeval distribution patterns of land and sea around the basins, the 'Xuefeng' paleo-uplift and the 'Qianzhong' paleo-uplift were the most possible provenance areas, this is because the Yangtze area in South China was confined by these ancient lands during this period [20]. In addition, the large-scale tectonic movement which led to the uplift of the Yangtze area [49] resulted in an overall regressive event in South China and the formation of a shallow-water deposition condition, under which the accumulation rate became accelerated. The increase in the arenaceous components within this interval also proves the shallowing of deposition water.

Redox Conditions
The redox condition of the sedimentary water controls the enrichment/depletion of certain elements, especially multivalent trace elements, and the preservation of total organic carbon (TOC) [50,51]. Thus, studying the redox conditions of the Wufeng and Longmaxi Fm. is critically important for shale gas exploration. Oxic conditions occurred Through the analysis above, it can be concluded that the provenance has changed several times during the depositional periods of Wufeng and Lower Longmaxi Fm. The Wufeng and the bottom of the Lower Longmaxi (83.5~101 m) are mainly characterized by biogenic inputs. During the deposition periods of middle and upper layers of the Lower Longmaxi, the strata are composed of terrigenous inputs (16~83.5 m). Until the topmost of the Lower Longmaxi (0~16 m), multiple provenance areas contributed deposition materials to the Basin simultaneously. The GYQ (Guanyinqiao) layer is a set of biogenic calcareous depositions.

Redox Conditions
The redox condition of the sedimentary water controls the enrichment/depletion of certain elements, especially multivalent trace elements, and the preservation of total organic carbon (TOC) [50,51]. Thus, studying the redox conditions of the Wufeng and Longmaxi The elements vanadium (V) and chromium (Cr) are typical multivalent trace elements [43,[53][54][55], which are highly sensitive to the redox conditions of bottom water. The ratio of V/Cr can reflect the redox condition of bottom water. Generally, V/Cr < 2, under oxic conditions and 2.0 < V/Cr < 4.25 in suboxic conditions, while under anoxic water, V/Cr > 4.25 [37].
On the Liutang section, the V/Cr ratio of all samples is between 1.05 and 9.28, with an average value of 2.42. As shown in Figure 4, the V/Cr ratio profile clearly indicates the redox conditions of the strata. The successions at 23~33 m, 62~73 m, 83.5~87 m, 91.5~93.5 m, and 95.5~97 m were deposited under suboxic conditions. The strata deposited in the intervals of 87~91.5m and 93.5~95.5 m were characterized by anoxic deposition conditions, while the remaining layers were deposited under oxic depositional conditions.
The Framboidal pyrite developed well in the whole section (Figure 6a-i). The morphological significance of the framboidal pyrite can also indicate the redox conditions of the underlying water during deposition. Generally, the framboidal pyrite formed under the oxic condition is called diagenetic pyrite, with characteristics of large and uneven particle sizes, while the framboidal pyrite formed under the reductive condition is called syngenetic pyrite with characteristics of small and even particle sizes. Wignall and Newton [56] pointed out that the maximum particle size of framboidal pyrite formed under the oxic-suboxic sedimentary water can be greater than 20 µm, and the maximum size of the framboidal pyrite formed under the reducing environment is generally less than 20 µm.

Chemozones
Chemostratigraphy can be defined as using the variability of the chemical elements, or fingerprints, to evaluate the stratigraphic relationships within a succession, and some essential proxies, including elements and ratios, can reflect the chemical variability [8,9]. As discussed above, the target successions show heterogeneities in terrigenous inputs, calcareous contents, redox, and paleoproductivity conditions. Thus, the proxies used in the present study for chemostratigraphic purposes are closely related to these aspects. The proxies used for the chemostratigraphic work and their geological meanings are tabulated in Table 2.  Figure 4 shows the chemostratigraphic scheme for the Wufeng and Longmaxi Fms. acquired based on the proxies tabulated in Table 2 and the constrained clustering analysis. For details, refer to chapter 'Chemostratigraphy zonation.' The division results of chemozones rely on the changing trends of the geochemical proxies.

Paleoenvironmental Significance of Chemozones
The proxies used for a chemostratigraphic scheme include biogenic/terrigenous inputs, redox conditions, paleoproductivity, and carbonate. The sedimentary paleoenvironment controls their contents. Thus, every chemozone implies a unique paleoenvironmental condition.
The chemozone C2 was almost all deposited in the shallow and oxic shelf environment, except C2-2A (dysoxic). At the bottom of C2 or C2-1A, the sequences turn from siliceous to aluminous, which indicates the first increase in terrigenous inputs for the whole section. This increase continues until C2-2A, which has typical characteristics with high values of SiO 2 /(SiO 2 + 5Al 2 O 3 ) (>0.5), V/Cr (2~4.25), TOC (>2%), and Si exc (>10). The terrigenous inputs rose again from C2-2B. This interruption in the increase in terrigenous inputs may be interpreted as the persistent biological influence of siliceous Radiolarian on C2-2A. Chemozones C2-3, C2-4, C2-5, and C2-6 are very similar geochemically, as all of them were deposited in a shallow shelf paleoenvironment. High calcium in C2-3 and C2-5 makes them distinguishable from other subchemozones (C2-4 and C2-6). Considering that the carbonate dilution function has not influenced the C2 (see Figure 5e), the calcium in C2-3 and C2-5 may originate from late diagenetic hydrothermal fluid rather than a biogenic source or seawater input, which is different from biogenic zones (C1-2 and C1-4).
Similar to C2, C3 is deposited in a shallow and oxic shelf environment. However, it is noted that the values of V/Cr and TOC are much lower in the latter, which means that C3 has a shallower and more oxidized sedimentary environment than C2. Furthermore, C3 is typically characterized by a fluctuating Si exc , SiO 2 /(SiO 2 + 5Al 2 O 3 ) (Figure 4). The terrigenous proxies, Rb/Al, Ti/Al, and Y/Al, show the same changing trends (Figure 7), indicating multiple provenances during C3 period.
The sedimentary paleoenvironmental revealed by various chemical proxies within chemozones can be attributed to global eustasy and paleoenvironments in South China. During the early and middle periods of the Upper Ordovician, the Yangtze area of South China belonged to the semi-restricted euxinic basin [23]. Wufeng shales, except for the top Guanyinqiao (GYQ) layer, were deposited under this anoxic and deep-water condition. Thus, C1-1 is typically characterized by anoxic, high TOC, and low calcium content. At the end of the late Ordovician period, global sea levels, including in South China, dropped dramatically under the effect of the Gondwana glaciation [22], resulting in the oxidization of seawater in the research area. Thence, C1-2 is characterized by more oxic and more calcium than C1-1. At the start of the Lower Silurian period, deglaciation occurred, and global sea levels rose due to global warming, which was partially ascribed to the volcanism [57]. Moreover, this event affected almost all of the worldwide epicontinental seas [58], including the Yangtze region of South China. Thus, C1-3 shows geochemical characteristics of anoxic, high TOC, low calcium, and even high paleoproductivity. Following this sea level rising event, the global sea level started to decline again [45], which contributed to regressive seawater and the formation of the prograding sequences in the Yangtze area of South China, and the research area turned into a shallow shelf [23]. The simultaneous global-warming paleoclimatic conditions [47,48] accelerated the weathering of surrounding paleo-uplifts leading to abundant terrigenous materials transported into the Longmaxi successions.
Consequently, subchemozones C1-4 and C2 are characterized by high terrigenous inputs, high oxygen, and low TOC. During the final deposition stage of the Lower Silurian Fm., the Yangtze area uplifted as a whole [49], and a complete marine regression occurred in the research area. Moreover, multiple paleo-uplifts, including 'Xuefeng' and 'Qianzhong,' brought massive terrestrial inputs into the basin, which helps explain why chemozone C3 has a multi-provenance characteristic. This extremely shallow-water condition also resulted in the lowest TOC and V/Cr in chemozone C3 of the section.

Regional Correlation
The global eustasy is the main factor controlling the deposition of the Wufeng and Longmaxi Fms., which is also validated by the geochemical proxies and previous research (e.g., [22,45,49]). Thus, the existing stratigraphic framework for the strata is acquired mainly based on the sequence stratigraphy. The wells JY-1, PY-1, and DY-1 (see Figure 1 for their detailed locations), located at the depositional center of Wufeng and Longmaxi, are actually the most-studied wells (e.g., [3,[59][60][61]), of which the JY-1 well is located at the Jiaoshiba gas field and is the first industrial gas field in China, with recoverable resources of approximately 200 billion cubic meters. Thus, we compared the chemostratigraphic scheme with the sequence schemes derived from these wells. Two kinds of stratigraphic schemes can be obtained from the sequence stratigraphic results for JY-1, PY-1, and DY-1 wells. The first is the results described by Wang et al. [60] that the Wufeng and the Lower Longmaxi are comprised of an alternating Transgressive System Tract (TST) and a Highstand System Tract (HST) (Figure 8). The second is described by Chen et al. [3], Guo [59], and Zhang et al. [61]: the Wufeng Fm. consists of TST and HST, while the Lower Longmaxi is composed of TST, Early Highstand System Tract (EHST), and Late Highstand System Tract (LHST) (Figure 8). The existing sequence stratigraphic schemes confirm that two major transgressive-regressive cycles have occurred during the Wufeng and Lower Longmaxi periods, which is consistent with our chemostratigraphic finding that the chemozones C1-1 and C1-3 represent two transgressive phases, and the other chemozones (C1-2, and C1-4 to C3) are the depositions during regressive phases.
The chemostratigraphic profile can be used to construct the sequence-stratigraphic framework (e.g., [7]). TSTs are defined by retrogradational successions with a constant increase in accommodation [52]. Thus, TSTs are characterized by a general decline in terrigenous proxies and basin restriction proxies (e.g., [7]). HSTs are defined by progradational successions with a constant decrease in accommodation [52] and are characterized by an increase in terrigenous inputs and minimal levels of basin restriction [7]. LSTs (Lowstand System Tracts) are also defined by progradational successions [52] but are characterized by an increase in proxies of basin restriction [7]. SiO 2 /(SiO 2 + 5Al 2 O 3 ) was used as a proxy of terrigenous inputs, lower SiO 2 /(SiO 2 + 5Al 2 O 3 ) indicates more terrigenous inputs. Accordingly, the chemozones C1-1 and C1-3 can be correlated to TSTs due to their high contents of SiO 2 /(SiO 2 + 5Al 2 O 3 ), V/Cr, TOC, Si exc , and low Ca, which indicate retrogradational successions. The chemozones C1-2, C1-4, and C3 can be correlated to HSTs due to their low contents of SiO 2 /(SiO 2 + 5Al 2 O 3 ), V/Cr, TOC, Si exc , and high Ca which indicate progradational successions. See Figure 8 for details.
Through the comparisons above, it can be concluded that chemostratigraphy is an efficient way of establishing stratigraphic schemes but with much higher resolution than traditional methods. For example, chemostratigraphy can recognize the boundaries of chemozones C1-4, C2-1A, C2-1B, C2-2A, C2-2B, C2-3, C2-4, C2-5, and C2-6, which is difficult for the sequence-stratigraphic method. Figure 8. The chemostratigraphic scheme and its comparison with the regional sequence-stratigraphic schemes. a is revised from Wang et al. [60]; b is from Guo et al. [59]; c is from Zhang et al. [61]; d is from Chen et al. [3].

Identification of Sweet Spot
The chemostratigraphic work above provides a precise stratigraphic scheme, and when combined with geochemical fingerprints, such as elements and their ratios, it can also help with the evaluation of stratigraphic properties [2,29,30], which is the preferred process for sweet-spot identification.
Generally, total organic carbon (TOC) content determines the designation of a layer as a sweet spot, and the content of brittle minerals (for example, quartz and carbonate) decides the possibility of shales being fractured. Based on the high-resolution division results above and the geochemistry analysis, it is clear that the chemozone C1 may be the optimal candidate for the sweet-spot layer due to its high TOC contents and anoxic-suboxic sedimentary conditions. Moreover, the abundant siliceous and calcareous matters (Figure 4), especially biogenic silica in C1-1, C1-3, and C1-4, and biogenic calcium in C1-2, deduce the difficulty in fracturing shales. Thus, chemozone C1 (83.5-101 m) is thought to be the sweet-spot layer for shale gas exploration and development.

Conclusions
The sedimentary geochemical analysis indicates that the macroscopically-homogeneous Wufeng and Lower Longmaxi shales show prominent geochemical heterogeneities Figure 8. The chemostratigraphic scheme and its comparison with the regional sequence-stratigraphic schemes. a is revised from Wang et al. [60]; b is from Guo et al. [59]; c is from Zhang et al. [61]; d is from Chen et al. [3].

Identification of Sweet Spot
The chemostratigraphic work above provides a precise stratigraphic scheme, and when combined with geochemical fingerprints, such as elements and their ratios, it can also help with the evaluation of stratigraphic properties [2,29,30], which is the preferred process for sweet-spot identification.
Generally, total organic carbon (TOC) content determines the designation of a layer as a sweet spot, and the content of brittle minerals (for example, quartz and carbonate) decides the possibility of shales being fractured. Based on the high-resolution division results above and the geochemistry analysis, it is clear that the chemozone C1 may be the optimal candidate for the sweet-spot layer due to its high TOC contents and anoxicsuboxic sedimentary conditions. Moreover, the abundant siliceous and calcareous matters (Figure 4), especially biogenic silica in C1-1, C1-3, and C1-4, and biogenic calcium in C1-2, deduce the difficulty in fracturing shales. Thus, chemozone C1 (83.5-101 m) is thought to be the sweet-spot layer for shale gas exploration and development.

Conclusions
The sedimentary geochemical analysis indicates that the macroscopically-homogeneous Wufeng and Lower Longmaxi shales show prominent geochemical heterogeneities in paleoproductivity, terrigenous inputs, and redox conditions. Generally, the upper Wufeng and the bottom of the Lower Longmaxi shales are characterized by high paleoproductivity, low terrestrial contribution, and anoxic deposition conditions, while the upper part of the Lower Longmaxi shales feature a gradual decrease in paleoproductivity, high terrigenous inputs, and oxic deposition conditions. The geochemical analysis also recognized two transgressive-regressive cycles on Liutang section, which is consistent with two coeval global eustasy events. Based on these geochemical heterogeneities, the chemostratigraphic scheme for the Wufeng and Lower Longmaxi shales was established using the constrained clustering analysis, and a total of 13 chemozones could be gained which are comparable to the regional sequence stratigraphic schemes and with higher resolution. Moreover, the chemozone C1 (83.5-101 m), characterized by high TOC content and high brittle minerals, is thought to be the most preferable layer for shale gas exploration and development.