Pore Structure Fractal Characterization and Permeability Simulation of Natural Gas Hydrate Reservoir Based on CT Images

Guangzhou Marine Geological Survey, China Geological Survey, Guangzhou 510075, China Southern Marine Science and Engineering Guangdong Laboratory, Guangzhou 511458, China Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China Center of Oil & Natural Gas Resource Exploration, China Geological Survey, Beijing 100083, China School of Energy Resources, China University of Geosciences, Beijing 100083, China Key Laboratory of Tectonics and Petroleum Resources, Ministry of Education, China University of Geosciences, Wuhan 430074, China


Introduction
Natural gas hydrate refers to one kind of ice-like crystalline compound formed by a methane-based hydrocarbon gas and water in a stable domain controlled by a certain temperature and pressure condition [1]. In the field of exploration and development of the natural gas hydrate, especially in hydrate classification, extraction, dissociation, experimental simulation, productivity evaluation method, etc., the work has become a research hotspot pursued by scientific research institutions, experts, and scholars [2][3][4][5][6][7][8]. The first offshore natural gas hydrate production test was conducted by China Geological Survey at Shenhu Area in South China Sea in 2017 [9]. During the depressurization process of production, the solid hydrate is decomposed into methane gas and water. The two-phase flow passes through the formation and then flows into the wellbore. In the initial stage, the flow rate is stable. However, as production progresses, the flow rate is going to be greatly reduced, resulting in a significant decrease in production. The above-mentioned multiphase flow within pore space is extremely complicated, and its seepage mechanism is not yet clear. It is speculated to be related to pore structure and mineral composition of the hydrate reservoir.
The differences of porosity, connectivity, and topology of pore space always determine the macroscopic seepage properties [10,11]. The pore space of the reservoir core can be obtained by multi-laboratory test methods [12][13][14][15]. The obtained pore space images and parameters can also lay the foundation for calculation of permeability [16][17][18]. Many numerical simulation methods [19][20][21] were used to obtain permeability of the reservoir, which greatly promotes understanding of the relationship between reservoir porosity and permeability, and the formulation of reservoir development programs.
Fractal geometry theory, as the powerful tool for quantification of complex space, has been used by many scholars in the study of characterization of the pore structure and modelling of transport phenomena in conventional and unconventional reservoirs. Yu and Cheng [22] established the fractal permeability model for the bidispersed porous media and introduced the algorithm of pore space fractal dimension and tortuosity fractal dimension based on the box-counting method. Cai et al. [23] established a fractal permeability model for creeping microstructures under different pressure conditions based on the data of water flooding experiment combined with CT scanning [24] for clayey silt sediments. Liu et al. [25] used fractal dimensions to establish a fractal theory-based relative permeability model of hydrate-bearing sediments. N'Diaye et al. [26] computed fractal dimension, lacunarity, and succolarity to provide a precise description of porosity and pore characteristics. Xia et al. [27] proposed that succolarity can be used to predict the permeability of low permeable sandstone reservoirs. Therefore, fractal geometry theory has been proved to be a good method to characterize the reservoir microstructure and to study the macroscopic differences in physical properties of reservoirs. In the study of clayey silt porous media, the macroscopic changes of reservoir physical properties can also be explained from the perspective of microstructure based on fractal geometry theory. However, few studies have used fractal quantitative parameters to characterize the microporous structure of clayey silt porous media and further used to study the change of permeability.
In this study, pore structure images of six hydrate reservoir samples were firstly acquired by micro-CT scanning. Then, the porosity, permeability, and fractal parameters were obtained based on the binarized pore images. Finally, how microporous structures affect the permeability of clayey silt reservoir samples was analyzed according to the relationship between permeability and pore structure parameters.

Materials and Facility
2.1. Hydrate Reservoir Samples. Six clayey silt samples are taken from the field of the first offshore natural gas hydrate production test, which is located at Shenhu Area in the South China Sea [9]. There are eight natural gas hydrate deposits, among which the W11-17 deposit (Figure 1) is selected as the optimal target for the production test. The lithology of the reservoir is clayey silt, of which the clay mineral content is high (26%-30%), mainly composed of montmorillonite and illite. The average median diameter of the sample particle is 12 μm. The average values of effective porosity and hydrate saturation are 33% and 31%, respectively. Pore space of the reservoir is filled with three phases: solid hydrate, free gas, and water. Detailed geological, ore body, seismic, and logging information of reservoir can be found in literature [9].

Methods
After the original CT scanning images obtained, a series of preprocessing including extraction, filtering, and binarization [28] is conducted to acquire the binarized pore space structure images. Based on the binarized images, a variety of data such as sample porosity, permeability, and fractal structural parameters can be acquired and then can be used

Geofluids
to construct a relational model between pore structure and physical properties.

Porosity Calculation.
Porosity ϕ refers to the ratio of the pore volume to the total volume of the rock. It is an important parameter to quantitatively describe the rock reservoir capacity. The expression is shown below: where V p is the pore volume, V b is the total volume of the rock, and V s is the total volume of the skeleton.
In microscopic sample images, the pore volume and the rock skeleton volume correspond to the pixel volume of different color regions. Since the size represented by each pixel is uniform in all images of the same resolution, Equation (1) can be transformed to where N p is the number of pore pixels, N b is the total number of pixels of the sample image, and N s is the number of sample skeleton pixels. The resulting binarized image is converted into an array matrix as displayed in Figure 2(b), and the sample image is represented by the pixel value of each pixel in the image. A sample image of a cube corresponds to a threedimensional array matrix. By counting the number "0" and "1" in the matrix, the number of pore pixels and rock skeleton pixels can be obtained, respectively.

Fractal Structural Parameters Calculation.
In the fractal geometry theory, fractal dimension, lacunarity, and succolarity are the main parameters of fractal analysis. Fractal dimension is used as the basic parameter of fractal quantitative characterization. Among plentiful research results, the boxcounting method is a widely used and accurate method by which the fractal dimension of pore distribution can be obtained based on the structural images of the tested sample.
In this study, all fractal dimensions are calculated based on the box-counting method as well. In the box-counting method, the image is covered by boxes with scale r, and the number of the boxes N covered with the pores is recorded. The expression is as follows: Therefore, the negative value of the slope of the linear fit for the box-number and scale of double logarithmic plots in the self-similar interval is the pore space fractal dimension D.
Lacunarity analysis is a popular analytical tool of fractal microstructure in recent years. It can quantitatively analyze the degree of fractal clustering and characterize the changes of natural images and can consequently distinguish different structures with the same fractal dimension. Lacunarity characterizes the distribution of the observed objects in the image and measures the frequency and magnitude of observed objects. The larger the lacunarity, the more intensive the distribution of the observed object.
The gliding box-counting method is a general measurement method of lacunarity. The value of lacunarity can be defined as a statistical moment function when q = 2 divided by the squared value of q = 1 [29,30]: Equation (5) is used to normalize lacunarity, which can exclude the influence of porosity.
Succolarity is one of the crucial parameters of fractal geometry, which indicates the ability of fluid flowing inside the medium, that is, the extent to which the fluid in the pore is permeable or flowing. It can be used to measure the connectivity of the patterns or structures in different directions.
The algorithm of succolarity is also based on the boxcounting method, which evaluates the fluid flow ability in different flow directions by applying a virtual pressure field. The value of succolarity can be characterized as the pressure coverage value divided by the maximum possible coverage value [31]: 3.3. Permeability Calculation. In order to study the relationship between complex microstructure and macroscopic seepage property of the gas hydrate reservoir, we selected the absolute permeability, which is associated with the sample structure, as the research object. Absolute permeability k of the rock core is the permeability measured when only one fluid (single phase) flows through the pores of the rock, the fluid does not react with the rock, and the flow of fluid is based on Darcy's law of linear percolation. When a singlephase fluid with viscosity μ flows through a porous medium whose cross-section is A, length is L, and pressure difference is ΔP, the amount of fluid Q flowing through the rock pores within a unit time is Based on Equation (7), the absolute permeability of the six samples was simulated using the finite volume method (FVM) [32]. Avizo® 9 software is employed to convert the six binary graphic models into FVM models for the absolute permeability simulation. Simulation work is done with the Avizo XLab Hydro Extension module which assumes that the fluid was laminar and the viscosity was 0.001 Pa·s. The velocity-pressure coupling solver is used, and the absolute convergence factor is 1 × 10 −4 . The fluid flow governing 3 Geofluids

Results and Discussion
(3) (1) (2) Figure 5: The pressure field (Pa) distributed on y direction of the permeability simulation of the six hydrate samples.

Geofluids
preprocessing procedures conducted, the binarized pore structure images of the six samples are obtained, which are expressed in Figure 3 (B1-B6).
In order to better illustrate the special physical properties of hydrate reservoir samples, we selected two conventional sandstone cores for comparison. The pixel number and pixel resolution for CT images of sandstone cores shown in Figure 4 are 512 3 and~1 μm. According to Figures 3 and 4, unlike the irregular intergranular pores of conventional sandstone reservoirs, the pore morphology of the hydrate reservoir mostly turns out to be spherical and uniform in the image. The reason is that the hydrate reservoirs are mainly composed of fine silt particles with good sorting, so that the intergranular pores are consequently regular. In addition, the hydrate reservoir is greatly affected by microbial processes during diagenesis, which leads to the existence of spheroidal bioclastics and interdental pores. Figures 3 and 4 obtained, porosity and permeability can be calculated based on these images. Figure 5 illustrates the pressure field distributed on y direction of permeability simulation of the six hydrate samples. Table 1 shows the porosity and absolute permeability of hydrate Samples 1-6 and conventional sandstone Cores a and b. The ranges of porosity and absolute permeability of eight samples and cores are 0.17~0.39 and 18.48~1497.73 mD, respectively. The simulation results on Core 1 x, Core 3 z, and Core 6 z express that the pressure field and the velocity field are opposite in direction, which results in the inability to calculate the permeability. And the cause of this condition is presumed to be the absence of connected pores on these directions.

Porosity and Permeability. After the binarized images shown in
According to the data in Table 1, the hydrate sample exhibits higher porosity and lower permeability compared with the conventional sandstone core, and the porosity and permeability are poorly correlated (as shown in Figure 6, the correlation coefficient is only 0.4669). That is, the porosity of Sample 1, 2, and 3 is 0.22, 0.26, and 0.17, respectively, but the absolute permeability is much smaller than that of sandstone Core a, whose porosity is 0.18. The porosity of hydrate Sample 4 and Sample 5 has reached a high level at about 0.3, but the permeability is only an average of 200 mD. When the porosity of conventional sandstone Core b is greater than 0.3, the permeability can reach more than 1000 mD. As for Sample 6, although the porosity is also approximately 0.2, its permeability is more than twice of that of Sample 1, 2, and 3. In summary, the permeability of the hydrate sample is lower than that of the conventional sandstone core and the distribution is complex, having poor correlation with the porosity. It is necessary to carry out further analysis for the microscopic pore structure of the sample by using fractal theory to find out the seepage mechanism. The range of 3D fractal dimension of eight samples and cores is 2.70~2.85. The fractal dimension characterizes a complex feature of the overall pore distribution, which is generally affected by many factors, especially porosity. The larger the porosity, the larger the fractal dimension, and the more complex the pore space distribution develops. From this perspective, in the case of low porosity, hydrate samples exhibit a high fractal dimension compared with conventional sandstone cores (Samples 1, 2, 3, and 6 > Core a and Sample 4 > Core b), which proves that the distribution of pore space of the hydrate sample is more complicated than that of the conventional sandstone cores (Table 2).

Lacunarity and Permeability.
In order to observe lacunarity (heterogeneity) of different samples and cores more intuitively, medians of the normalized values at different scales are listed in Table 3. The range of lacunarity midvalue of eight samples and cores is 0.168~0.447. Among the hydrate samples with smaller porosity, it can be seen that though the porosities of Sample 1, 2, 3, and 6 are similar, permeability of Sample 6 is almost double of that of Sample 1, 2, and 3. The results indicate that with the heterogeneity of sample increases, the permeability decreases. The reason is that when porosity is small, permeability of the hydrate  Figure 6: Fitting schematic diagram of porosity and permeability of six hydrate samples (permeability is taken on y direction). 6 Geofluids sample with uniform pore distribution depends on the connectivity of pore throats. That is, the stronger the pore distribution heterogeneity of the hydrate sample, the worse the connectivity of the pore throats, and the smaller the permeability in consequence. Therefore, for a hydrate sample with a small porosity, the smaller the calculated lacunarity (the weaker the heterogeneity), the easier it is to exhibit a larger permeability. That is, the stronger the pore distribution heterogeneity of the hydrate sample, the worse the pore throat connectivity, and the smaller the permeability consequently. Therefore, for a hydrate sample with small porosity, the smaller the calculated lacunarity (the weaker the heterogeneity), the more likely it is to possess higher permeability. On the other hand, hydrate samples with larger porosity and sandstone cores show that the stronger the heterogeneity, the higher the permeability. Namely, Core a is more heterogeneous in the case of possessing smaller porosity, and its permeability is basically not lower than Samples 4 and 5; Sample 5 has a lower porosity than Sample 4, but it is highly heterogeneous, and the permeabilities in two directions are greater than Sample 4. The reason is that in the hydrate samples with larger porosity and sandstone cores, the distribution of pore space has already formed pore throats for fluid flow. Currently, the greater the heterogeneity of the pore distribution, the larger the diameter of the pore throats, which consequently leads to the higher permeability. Therefore, in hydrate samples with large porosity and all sandstone cores, the larger the calculated lacunarity (the stronger the heterogeneity), the greater the permeability. Table 4 shows the calculated value of the 3D succolarity on six directions of six hydrate samples. The range of the results is 0.0024~0.3312. It can be seen that, except for Sample 3, hydrate samples have little difference in succolarity values on all directions. Since succolarity can characterize the physical significance of the anisotropic characteristics of core pore structure distribution, the hydrate sample pore structure has no obvious anisotropic characteristics. For Sample 3, the extremely small succolarity value on z direction also indicates that there is no connected pores in this direction, which explains why the permeability cannot be calculated. Figure 7 is a schematic diagram showing the fitting of succolarity and permeability on different positive directions of six hydrate samples (except Sample 1 x, Sample 3 z, and Sample 6 z, on which direction the permeability cannot be calculated). It can be seen that, succolarity values of the hydrate samples show a good correlation with the permeabilities (correlation coefficient is 0.663). The cause of the incomplete fitting is that the pore radius of Sample 5 is larger ( Figure 5), which results in a higher permeability, however    7 Geofluids the effect of pore radius has not been considered by succolarity. Even so, compared with porosity, succolarity can be used to fit the permeability and establish a prediction model, which gives a certain theoretical guidance to the changing rule of hydrate reservoir permeability.

Conclusions
In this study, six natural gas hydrate samples are conducted with CT scanning to get the pore structure, which are used to calculate the porosity, permeability, and fractal parameters. The results show that the fractal parameter can quantitatively characterize the pore structure and analyze the change of permeability of the hydrate reservoir. The following conclusions can be derived: (a) Compared with the conventional sandstone cores, the hydrate samples show higher porosity and lower permeability, and the correlation between porosity and permeability is poor (b) From the comparison results of the fractal dimension, hydrate samples possess higher fractal dimension than conventional sandstone cores and the structural development is more complicated (c) From the results of lacunarity calculation, it turns out that for hydrate samples with smaller porosity, with the lacunarity of pore distribution (the weaker the heterogeneity) decreases, the permeability increases; for the hydrate samples with larger porosity and sandstone cores, with the lacunarity of pore distribution (the stronger the heterogeneity) increases, the permeability increases (d) According to calculation results of succolarity, hydrate samples basically show nonanisotropic development characteristics. It is clear that the correlation between the succolarity value and the permeability value is high. Therefore, a predictive model of permeability can be built to provide theoretical support for hydrate reservoir development

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.