The Description of Shale Reservoir Pore Structure Based on Method of Moments Estimation

Shale has been considered as good gas reservoir due to its abundant interior nanoscale pores. Thus, the study of the pore structure of shale is of great significance for the evaluation and development of shale oil and gas. To date, the most widely used approaches for studying the shale pore structure include image analysis, radiation and fluid invasion methods. The detailed pore structures can be studied intuitively by image analysis and radiation methods, but the results obtained are quite sensitive to sample preparation, equipment performance and experimental operation. In contrast, the fluid invasion method can be used to obtain information on pore size distribution and pore structure, but the relative simple parameters derived cannot be used to evaluate the pore structure of shale comprehensively and quantitatively. To characterize the nanoscale pore structure of shale reservoir more effectively and expand the current research techniques, we proposed a new method based on gas adsorption experimental data and the method of moments to describe the pore structure parameters of shale reservoir. Combined with the geological mixture empirical distribution and the method of moments estimation principle, the new method calculates the characteristic parameters of shale, including the mean pore size (x¯), standard deviation (σ), skewness (Sk) and variation coefficient (c). These values are found by reconstructing the grouping intervals of observation values and optimizing algorithms for eigenvalues. This approach assures a more effective description of the characteristics of nanoscale pore structures. Finally, the new method has been applied to analyze the Yanchang shale in the Ordos Basin (China) and Longmaxi shale from the Sichuan Basin (China). The results obtained well reveal the pore characteristics of shale, indicating the feasibility of this new method in the study of the pore structure of shale reservoir.


Introduction
Shale has attracted increasing attention recently due to its nanoscale pores containing abundant oil and gas resources and has become a focus of global unconventional oil and gas exploration and development. Compared to conventional reservoirs, shale is a tight reservoir with a complex origin and strong heterogeneity [1][2]. Their pores have a complex geometry at the nanoscale and have low porosity and ultra-low permeability [3][4][5][6]. Previous studies have investigated the pore size and distribution of different shale formations. The Barnett shale in North America has various types of pores [6][7][8], with an organic nanopore distribution between 5 and 750 nm, a median size of approximately 100 nm [6]. The pore size of the marine shale in South China is between 2 and 900 nm, with an average of approximately 3 to 40 nm [9][10][11][12][13][14]. The pores of the continental shale of North China have a distribution between 2 and 35 nm [15][16][17][18]. The abundant nanoscale pores of shale provide good spaces for gas accumulation, playing an important role for gas storage and migration [1,6,8,[19][20][21]. Thus, the study of pore structure features is of great significance for shale oil and gas evaluation as well as for exploration and development.
Pore structure refers to the geometry, size, distribution and the interconnected relationships between the pores and throats in rocks [22]. In recent years, transmission/scanning electron microscopy (TEM/SEM), field emission scanning electron microscopy (FE-SEM), atomic-force microscope (AFM), nuclear magnetic resonance (NMR), small-angle and ultra-small-angle neutron scattering (SANS and USANS), X-ray micro tomography (XMT), mercury injection capillary pressure (MICP) and gas adsorption have been used to study the pores of shale [3][4][5][6][7][8]20,[23][24][25][26][27][28]. Image analysis methods can be used to observe the pores developed in rocks through TEM/SEM, FE-SEM and other high resolution microscopies. The obtained image data can be used to analyze the geometry and connectivity of pores. Although a large quantity of pore throat information can be acquired intuitively using image analysis, the results obtained depend on the device resolution, and a large number of experiments are needed due to the small imaging scope. Moreover, it is difficult to conduct quantitative analyses with image analysis. NMR, SANS, Nano-CT and other radiation methods be used to conduct non-destructive analysis and acquire detailed data rapidly, but the measurement errors are strongly influenced by sample preparation and heterogeneity. In addition, radioactive sources are rare and quite expensive. Researchers often use a fluid invasion method (e.g., MICP, N 2 adsorption, or CO 2 adsorption) to determine the shale pore size distribution and use theoretical models to calculate the specific surface area, the volume and other parameters.
In the gas adsorption method, the volume of liquid nitrogen filled in test samples is an equivalent volume of the pore space under isothermal conditions. A relationship between relative pressure (P/P 0 ) and gas absorption (desorption) volume (V) can be constructed by measuring gas adsorption volumes under different relative pressures (P/P 0 ). The resulting adsorption and desorption isotherms can be used to obtain the specific surface area and pore size distribution using theoretical model calculations [29]. At present, commonly used theoretical models include Barrett-Joyner-Halenda (BJH) [30], Horváth-Kawazoe (HK) [31], density functional theory (DFT) [32] and quenched solid density functional theory (QSDFT) [33][34] models. According to the pore size classification standard of the International Union of Pure and Applied Chemistry (IUPAC), the three categories of pores are micropores (diameter < 2 nm), mesopores (2 nm < diameter < 50 nm) and macropores (diameter > 50 nm) [35]. Of the models mentioned above, the BJH model can effectively measure the mesopores with pore sizes of 2 to 50 nm [30,36]. The HK model is commonly used to analyze microporous materials [31]. The DFT model is based on the Tarazona equation and QSDFT model, which improves upon the DFT model and can be used to effectively analyze microporous-mesoporous materials [32][33][34]. Although it has been extensively used to characterize the pore structure of shale, the gas adsorption method still has some drawbacks: e.g., it requires complex theoretical model calculations, and the results obtained by different models show larger errors [37][38].
Theoretically, the gas adsorption method can measure open pores with a diameter larger than 0.35 nm, which is the diameter of nitrogen molecules. However, only pores with a diameter between 1.5 and 200 nm can be measured effectively due to the limitations of instrument accuracy and experimental operation difficulty. Furthermore, the parameters obtained by the gas adsorption method are relative simple and cannot be used to evaluate the mean pore size, pore sorting and pore structure comprehensively and quantitatively.
The major goal of this paper is to proposed a new method to describe the pore structure parameters of shale reservoir. The new method is based on gas adsorption experimental data and, combined with the geological mixture empirical distribution and the method of moments estimation principle, optimizes eigenvalue algorithms and determines pore structure parameters by applying the method of moments. Using this method, the calculation of the parameters of the pore structure is simplified, and the nanopores in shale can be described. This study aims to expand the current methods for evaluating the pore structure of shale, and simplify the algorithm of eigenvalues.

Method Description
Based on extensive investigations Luo and Wang (1981) discovered that the pore throat distribution of limestone and carbonate rocks did not comply with the normal distribution proposed by Chilingar (1972). Instead, the distribution was a combination of various pore throat distributions, which should be caused by multiple factors of diagenesis and epidiagenesis. For a mixed distribution, the eigenvalues can be determined by the method of moments according to the numerical characteristics of the empirical geological distribution [22]. Using the method of moments, Luo (1981) and Wang (2003) extensively studied the low permeability reservoirs of the oil and gas bearing basins in China, such as the Ordos Basin, Sichuan Basin, and Bohai Bay Basin. The characteristic parameters obtained revealed the pore structure characteristics of reservoir rocks and provided constructive guidance for the exploration and development of oil and gas fields [22,39].
The obtained date should be grouped before the method of moments is used to calculate the characteristic parameters of pore structure. The shale reservoirs with low permeability, low porosity of nanoscale pores belong to the tight-ultra-tight reservoirs [3][4][5][6]. Thus, the grouping interval of the observation values of sandstone and carbonate rocks is not suitable for the grouping of shale observation values, and a new grouping observation interval is required. By taking the gas adsorption volumes (under liquid conditions) during adsorption experiments as the observation values (equivalent to saturation), the functional relationship between relative pressure (P/P 0 ) and mean pore size (F) can be established by the classical Kelvin equation [40] Halsey equation [29] t ¼ 0:354 ½À5 = lnðP=P 0 Þ 1=3 ð2Þ and conversion equation [22] F where r k is the capillary radius in nm, P/P 0 is the relative pressure, t is the multilayer thickness in nm, r p is the pore radius in nm and D is the pore size (diameter) in nm in Eq (3) and expressed by mm in Eq (4). According to the Oil and Gas Industry Standard of People's Republic of China [41] and the National Standard of the People's Republic of China [42] for gas adsorption, the relative pressure range (P/P 0 = 0.250-0.995) determines the scope of F value: 11.35-18.35F. With equal spacing of F, the observation values can be divided into fifteen intervals with a width setting of 0.5F. This classification method actually shortens the coarse pore part intervals, widens the fine pore part intervals and considers the proportion of different magnitude pore weights. As such, the calculated characteristic parameters reveal the pore characteristics of shale well. The method of moments can analyze and evaluate the pores more effectively with the refined interval.
The observation values can be summarized as a mathematical process using mathematical language. For many types of observation values, important mathematical characteristic parameters of rock pore throats include: The mean is one of the parameters to describe the centralized location of data; it describes the average position of experimental data. It represents the average value of the entire pore size distribution for the pore structure of reservoir rocks. This value can be acquired through the weighted average of the observation values in our study, namely where x i is the start value (mid-value or end-value) of the i-th interval, which can be expressed by F for reservoir rocks. A greater F value indicates a smaller pore size (D). f i is the observation value (ΔS i ), i.e., the nitrogen saturation under liquid conditions that can be expressed as a percentage, S max is the maximum inlet nitrogen saturation and n is the number of groups with a default value of 15. This algorithm accounts for the proportion of pores of different sizes according to the weighted data and is more in line with the realistic situation of pore size distribution.

Standard deviation (σ)
The standard deviation is the scattered characteristic parameter. It is the mean-centered measure of the degree of dispersion of the data. It describes the degree of dispersion of the experimental data over the entire number axis and represents the extent of dispersion from the mean value. In this paper, the standard deviation refers to the degree of dispersion of pore size with respect to the mean and can also be called the sorting coefficient of the pores. For the pore system, a smaller standard deviation of pore indicates a better sorting coefficient.
3. Variation coefficient (c) The variation coefficient, defined as the ratio of the standard deviation to the mean value, is a useful measurement for the variation of observation values. It can be applied to describe the ratio between the mean pore size (F value) and sorting degree. If the mean pore size (F value) is large (more fine pores) and the sorting degree is good (a smaller degree of dispersion for the pores), the value of c will be small. Within a certain range, c values can be used to evaluate the quality of the reservoir pore structure. In general, a large value of c indicates a more discrete degree of the pores. A large difference in the pores size is conducive to gas storage and migration, indicating a good pore structure of the reservoir rock.
The skewness is a measurement of the direction and extent of the data distribution. In geology, skewness reflects the pore size distribution; it describes the inclination to a larger or smaller size from mean value, with a coarse skewness indicating larger pores, and vice versa.
In [22] Luo treated the S max in Eqs (5), (6) and (8) as a constant of 100 by regarding the maximum mercury saturation as 100% during the MICP. Li (2001) suggested the calculated results would be overestimated when taking the S max as a constant 100, which is impractical [43]. Noticed that shale has low porosity and ultra-low permeability and that closed pores may exist within shale. Therefore, pores may not always be fully filled by nitrogen, and saturation may not reach 100%. As such, in our method, S max was introduced into the denominator of Eqs (5), (6) and (8)

Samples and experiments
Twenty-two samples (Y-1 to Y-22) were collected from the Triassic Yanchang shale in the southern area of the Ordos Basin (China), and fifteen samples (Y-23 to Y-37) were collected from the Silurian Longmaxi shale in the southeast area of the Sichuan Basin (China). Nitrogen adsorption was conducted in the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation of Chengdu University of Technology, Chengdu, China. The QUADRASORB SI specific surface and porosity analyzer was used to test the thirty-seven shale samples and obtain the nitrogen adsorption data.
The relative pressure range of the instrument is 0.005-0.995, and the measuring range is 0.35-400 nm. The samples crushed to less than 0.250 mm. Before analysis, every sample should be degassed by vacuum to ensure that adsorbed moisture and volatile matter were removed. In our study, samples were dried under a vacuum at 80°C for at least 22 h. After drying, samples underwent the nitrogen adsorption test under the relative pressure range from 0.005-0.995 at 77.3K.

Data processing
Sample Y-1 illustrates the data processing process. First, an adsorption-desorption curve was plotted (Fig 1A) from the original data from the nitrogen adsorption experiment ( Table 1). The nitrogen desorption volumes V i (g) under different relative pressures were extracted from the desorption branch and then converted into liquid volumes V i (l) to calculate the inlet nitrogen saturation (S i ) and interval saturation (ΔS i ). Finally, the characteristic parameters, such as mean pore size (x), standard deviation (σ), variation coefficient (c) and skewness (S k ) can be calculated by substituting S i and ΔS i into Eqs (5)(6)(7)(8), as shown in Table 2. The experimental data of the other thirty-six shale samples were processed by the above method, and the results are shown in Table 3.

Results
N 2 isotherms demonstrate a wide range in adsorption and show a hysteresis pattern for all samples (Fig 1). The shape of the hysteresis loop reflects the pore geometry [20,35,44]. The adsorption isotherms show an upward trend in areas of high pressure, implying the existence of mesopores and macropores [35,45]. In the medium-pressure area, the adsorption curve does not coincide with the desorption curve, forming a hysteresis loop, indicating a capillary  3  condensation phenomenon [35,40]. A capillary condensation phenomenon indicates the existence of open pores, e.g., ink bottle-shaped pores, slit-shaped pores and conical pores [35,44]. According to IUPAC [35], the hysteresis loop of thirty-seven samples mainly falls into two types. Samples Y-1 to Y-13 and Y-23 to Y-31 can be classified as type H2 (Fig 2; ink bottleshaped pores), and samples Y-14 to Y-22 and Y-32 to Y-37 can be classified as type H3 (Fig 2; slit-shaped pores). Sample ID The pores in the shale of the Yanchang Formation have a mean pore size distribution between 14.66 and 17.45F, i.e., 5.59-38.63 nm, a standard deviation ranging from1.40 to 2.00, and a variation coefficient exhibiting a narrow distribution from 0.081 to 0.126. The distribution of skewness is from -1.930 to 0.302, dominated by negative contributions, implying the pore size (D) is smaller than the mean pore size (D) (#Y-1 to #Y-22 in Table 3). The mean pore size of Longmaxi shale is between 16.45 and 17.98F, i.e., 3.87-11.17 nm, with a standard deviation from 0.76 to 2.17, a variation coefficient from 0.042 to 0.175, and a skewness distribution of -3.444 to -0.698 (#Y-23 to #Y-37 in Table 3).
In accordance with the IUPAC classification [35], the pores from Yanchang shale and Longmaxi shale are mesopores. According to the principles of geological mixture empirical distribution and method of moments, the standard deviation reflects pores sorting, the variation coefficient reflects whether the pore structure is good or not, and the skewness describes whether the pore size tends toward coarse or fine [22,39]. The parameters derived from the method of moments directly reveals that the shale mean pore size (diameter) is small, the pores well sorted and the pore structure is poor. These parameters indicate poor storage and migration conditions.

Discussion
According to the calculated characteristic parameters in Table 3, the relationship between the mean pore size (F) and standard deviation (sorting coefficient), variation coefficient and skewness were established (Fig 3). As seen from Fig 3A, the mean pore size (F) and standard deviation of Yanchang shale has a certain correlation (R = 0.467, R 2 = 0.218): a smaller F value is associated with a larger pore size (D) and poorer pores sorting, it reflects the pores of larger pore size (D) has a large proportion of the total pores. In contrast, a larger F value, is associated with a smaller size (D) and better pores sorting, i.e., a large proportion of pores of smaller pore size (D). The mean pore size (F) of Yanchang shale shows a significant relation to the variation coefficient (R = 0.730, R 2 = 0.533), as shown in Fig 3B. A small F value, corresponds to a large pore size (D), which is favorable for oil and gas storage and migration, suggesting a good pore structure (large variation coefficient). In contrast, if the F value is large, the pore size (D) is small, which is not favorable for oil and gas storage and migration, suggesting a poor pore structure (small variation coefficient). The skewness shown in Fig 3C extends from -1.930 to 0.302, which is dominated by the negative contribution, indicating that the pore size (D) is smaller than the mean value (D), i.e., the pore size distribution of most samples is inclined to fine skewness.
Similar to Yanchang shale, the standard deviation and variation coefficient of Longmaxi shale also showed a significant relationship with the mean pore size (F) (R 2 = 0.839 in Fig 3D, R 2 = 0.900 in Fig 3E). With increases in the F value, the standard deviation decreases and the same variation coefficient decreases, and vice versa. The pore size distribution of most samples is inclined to a value smaller than the mean value (D) (Fig 3F).
Compared with sandstone and carbonate reservoirs, shale reservoirs are tighter and, having experienced a series of diagenesis, have a stronger heterogeneity and possess more complex pore structures, with a nanoscale pore system being dominant [1,[3][4][5][6]44]. According to the numerical characteristics of the geological mixture empirical distribution, the method of moments, which accounts for the diagenesis of reservoir rocks and the impact of epidiagenesis on rock pore structure [22], can be used to determine the characteristic parameters of the shale pore structure. To discuss the reasonableness of the calculated results and explore the feasibility of studying the shale pore structure using the method of moments, the calculation results have been compared with those calculated by gas adsorption theory. N 2 adsorption isotherms indicate a wide adsorption range and show a hysteresis loop for all samples (Fig 1). Samples Y-1-Y-13 and Y-23-Y-31 can be classified as type H2, and samples Y-14-Y-22 and Y-32-Y-37 showed type H3 (Fig 2).  stated that pores of type H2 are favorable for gas storage but not for migration, whereas the pores of type H3, having both well-developed micropores and macropores, are favorable for gas migration due to the good pores connection [4,45]. Following this viewpoint, the standard deviation and variation coefficient of samples with type H2 are smaller than type H3. As seen from Table 3, the results of the method of moments analysis show that samples Y-1 to Y-13 and samples Y-14 to Y-22 have pore sizes (D) of 5.59 to 13.75 nm and 18.79 to 38.63 nm, respectively. The standard deviations for the two types of samples are 1.40 to 2.00 and 1.60 to 1.96, respectively. The former have a small pore size (D) and are well-sorted. The variation coefficients of 0.081 to 0.122 and 0.096 to 0.126 for the two types of samples indicate the latter samples generally have a larger value and are better for gas storage and migration. Finally, the results show the skewness of -1.930 to 0.035 and -0.076 to 0.302 for the two types of samples. Whereas samples Y-1 to Y-13 have a smaller pore size (D), better sorting, and a smaller variation, they are poorer for gas migration compared with Y-14 to Y-22. The characteristic parameters of the Longmaxi shale sample also indicate that the type H3 hysteresis loop of rock samples are conducive to gas storage and migration compared with type H2 samples (Table 3). Thus, the analysis results of the method of moments are consistent with the pore characteristics revealed by isotherm hysteresis loops, suggesting that the characteristic parameters derived from the method of moments can effectively reflect pore structure characteristics.
The mean pore size (D) derived from the method of moments is discrepant with that derived from the BJH method (Table 3). This discrepancy, may be due to the following reasons: The application of the BJH method depends on two fundamental assumptions: (1) the pores are cylindrical in shape, and (2) the physical adsorption occurs on the pore walls and capillary condensation occurs in the inner capillary volume [30]. The pore structure of shale is highly complex, with various types of pores, such as cylindrical pores, narrow slit-like pores, conical pores and ink bottle-shaped pores [20,35,45]. Shale does not have only cylindrical pores. Thus, the application of the BJH method might result in relatively large errors in the estimated pore size. Additionally, the primary purpose of the BJH method is to calculate the pore size distribution. The average pore diameter is estimated by assuming a cylindrical model [46] without considering the pores weighting within different pore size ranges. It does not conform to the actual situation in shale and makes the results less meaningful. In contrast, considering the characteristics of pore sizes from the entire range, the mean pore size calculated by the method of moments is derived from the weighted observation values and can comprehensively reflect actual shale reservoirs.
Summarizing the analysis and discussion above, the characteristic parameters calculated by the method of moments are consistent with the actual characteristics of the shale pore structure and can describe the pore structures of shale reservoir in more detail. This result suggests that studying the pore structure of shale by the method of moments is feasible and that the result is reasonable.

Conclusions
(1) The method of moments considers the influence of diagenesis and epidiagenesis on the shale pore structure, which is consistent with the geological regularity of pore size distribution and can be used to describe pore structure characteristics through various parameters quantitatively and accurately. Thus, the method of moments can be used to study the pore structure of shale.
(2) The grouping interval of shale observation values was reconstructed based on the established functional relationship between the relative pressure (P/P 0 ) and mean pore size (F). By determining the observation value interval to be 11.35-18.35F and taking 0.5F as the spacing, fifteen intervals were created, and the algorithm of eigenvalues was optimized. The case study indicated that the pore structure characteristic parameters calculated by the method of moments were reasonable and reliable when compared with those derived from gas adsorption theory. This result suggests that the new method could be used to quantitatively describe pore structure characteristics.
(3) The mean pore size (F) of the shale reservoirs of the Yanchang Formation in the Ordos Basin was related to the standard deviation: a larger F value indicated better sorting, and vice versa. The mean pore size (F) was significantly related to the variation coefficient: a larger F value indicated a smaller variation coefficient and poorer pore structure, and vice versa. Similar to Yanchang shale, Longmaxi shale also exhibited these relationships.
(4) The mean pore size of Yanchang shale is 14.66-17.45F, the standard deviation distribution is 1.40-2.00, the variation coefficient ranges from 0.081-0.126 and the skewness is -1.930 to 0.302. The mean pore size of Longmaxi shale is 16.45-17.98F, the standard deviation distribution is 0.76-2.17, the variation coefficient ranges from 0.042-0.175 and the skewness is -3.444 to -0.698. Yanchang shale and Longmaxi shale both have a small mean pore size (diameter), nice pore sorting, a small pore size (diameter) degree of dispersion and fine skewness.