Influence of Particle Size Distribution on the Physical Characteristics of Pore-Filling Hydrate-Bearing Sediment

Nature gas hydrates (NGHs) are regarded as a potential alternative energy source due to their huge reserves and wide distribution. According to the geophysical surveys, the pore-filling hydrates occupy a large proportion of the global hydrate reserves, especially for the marine regions. Therefore, with a novel pore-scale 3D morphological modeling algorithm, this study systematically studied the effect of the particle size on the physical characteristics of the pore-filling hydrate-bearing sediment (HBS). The pore system evaluations and permeability simulations were performed by utilizing pore network modeling (PNM), and the thermal and electrical simulations were conducted by utilizing a finite volume method (FVM). The results show that for the HBS with smaller particle size, the average radius of the pores and throats would also be reduced, and the fractal dimension of the pore system would be increased. In addition, with the increasing hydrate saturation, the fractal dimension of the pore system will increase firstly and then decrease. And these parameter evolutions could impact the physical properties correspondingly; specifically, the decreasing particle size in the HBS would reduce the permeability and electrical conductivity of HBS and enhance the apparent thermal conductivity of HBS.


Introduction
Nature gas hydrates (NGHs) are crystalline inclusion compounds which can be formed under the condition of low temperature and high pressure [1]. Due to the huge reserves and high energy density of the NGH, it is considered as an alternative energy source [2]. There are mainly three occurrence forms of the NGH: (1) the pore-filling hydrate which exists in the pore space, (2) the load-bearing hydrate which contacts the adjacent mineral particles, and (3) the cementing hydrate which exists between the mineral particles [3]. However, according to the recent results of the geophysical surveys, the pore-filling hydrates could be directly observed or deduced. This is because the ocean is mainly excess-water environment, and under the water-containing condition, hydrates are prone to nucleate in the middle of the pores to form pore-filling hydrates [4]. Therefore, pore-filling hydrates occupy a large proportion of the global hydrates in deep sea [5].
The physical characteristics of the hydrate-bearing sediment (HBS) are important to develop a reasonable exploitation plan [6,7], which mainly include the permeability [8], the thermal conductivity [9], and the electrical conductivity. In terms of HBSs' permeability research, Li et al. [10] measured the permeability of unconsolidated quartz sands with different particle sizes before and after hydrate formation and found that the larger the particle sizes of the quartz sands, the higher the absolute permeability will be. Kumar et al. [11] measured the glass-bead porous medium's absolute permeability, and the permeability attenuation index was revised according to the measured data. In order to further study the influencing factors of permeability, Yang et al. [12] combined the X-ray computed tomography (CT) to study the effect of different pore and throat radius, hydrates' distribution, and fluid flow direction on permeability of the HBS. Wang et al. [13] utilized CT to study the porous medium anisotropy's influence on the relative permeability of water and gas during the formation of hydrates and found that the relative gas permeability was higher in the horizontal direction, while the relative water permeability was higher in the vertical direction. In terms of HBSs' thermal conductivity research, the thermal conductivity of the NGH in the permafrost region was measured and was compared with that in the ice, finding that in the NGH was lower than that in the ice [14][15][16][17][18][19]. Some researchers utilized the single-sided transient plane source technique to measure the thermal conductivity of HBS at a certain temperature and pressure and then proved the reliability of the molecular dynamics simulation according to these data [20]. In addition, Muraoka et al. [21] utilized the hot-disk transient method to measure the thermal conductivity of the NGH and that of the mud layer samples from the Nankai Trough. In terms of HBSs' electrical conductivity research, some researchers measured the electrical resistivity of the glass bead specimen [22] and found that when the saturation index was constant, the Archie equation could not describe the relationship between the hydrate saturation and electrical resistivity. Moreover, Liu et al. [23] simulated the formation of hydrates in the quartz sand and found that the electrical resistivity would increase with hydrate saturation increasing.
However, the previous laboratory tests cannot strictly synthesize the HBS with the same occurrence form to study its physical characteristics. At the same time, due to various geological and historical reasons, the particle size of the reservoir has a wide distribution range. For this reason, in this study, 24 pore-filling HBSs with different saturations and particle sizes are generated by combining a real digital core and a hydrate modeling algorithm. And then, these factors' effects on the evolution of the pore characteristic, the perme-ability, the thermal conductivity, and the electrical conductivity of the HBS are analyzed.

Materials and Facility.
In this study, we utilized the three different particle sizes of Fujian standard sands (ISO Co., China) to generate the host sediment and named them ST1, ST08, and ST06 separately. The particle size distribution diagram of these sands is shown in Figure 1. Then, these sands were placed in a measuring cylinder whose diameter was 20 mm and the tamper was utilized to compress so that the specimen could be shaped. This study utilized an X-ray CT (SMX-225CTX-SV, Shimadzu Co., Japan) device to perform the imaging of the host sediment. The relevant parameters about this device have been described in detail by Yang et al. [12]. After 360°rotations of a specimen, a series of two-dimensional slice images were collected and a median filter was utilized to process these two-dimensional images. Then, we utilized the inspeXio (Shimadzu Co., Ltd, Japan) to reconstruct these two-dimensional CT images into a three-dimensional model and chose the middle area as the representative elemental volume (REV) according to the operation of Wu et al. [24]. The REV is a cube consisting of 500 × 500 × 500 voxels, and the voxel size is 21 μm/voxel.

Methods.
In recent decades, many modeling algorithms have been proposed for reconstructing the digital core of different sediments, such as multiple-point statistics [25] and the genetic algorithm [26]. However, they cannot be used  for remold hydrate in the digital core, let alone controlling the occurrence type saturation accurately.
A novel pore-scale 3D morphological modeling algorithm for hydrate generation proposed by Wu et al. [24] was used in this study. As illustrated in Figure 2(b), the spe-cific steps are as follows: Firstly, conducting CT scanning for the host sediment. Secondly, segmenting the sand and pore phases according to the gray value difference. Thirdly, conducting the erosion operation to the image stacks of pore phase with certain voxels. Fourthly, conducting the dilation 3 Geofluids operation to the eroded image stacks of pore phase with certain voxels, and the obtaining image stacks is the image stacks of the pore-filling hydrate. Finally, assembling the sand phase, hydrate phase, and corresponding pore phase for the image stacks of the pore-filling HBS. The more detailed introduction could be found in the morphological algorithm proposed by Wu et al. [24]. In addition, the bulk variation since the hydrate generation is not considered, this is because Lee et al. [27] found that during hydrate generation, the HBS's volumetric strain is only around 0.1% which could be insignificant. After these operations, 24 models were generated, as can be seen in Figure 2. In these models, the yellow area represents the sand, while the gray area represents the hydrate. In addition, it could be found that the radius of hydrate particles in pore space would increase gradually as the hydrate saturation increases, and some hydrate particles that do not contact with the sand particles (pore-filling type) will gradually contact with the sand particles (load-bearing type), which is consistent with previous studies [28,29]. Specifically, hydrate saturation calculation was as follows: where S h is the hydrate saturation, V h is the volume of hydrate, and V ps is the volume of the pore space without hydrate. Then, the pore network modeling (PNM) and finite volume method (FVM) were utilized to simulate the effect of the particle sizes and hydrate saturations on the HBS's physical characteristics.

Results and Discussion
3.1. Characteristic Evolution of Pore Space. The pore space could impact the reservoir capacity significantly [29,30], and the PNM is an effective tool to evaluate pore space evolution. In this study, PNM was utilized to estimate the pore and throat characteristics in the porous media under different particle sizes. Figure 3 shows that as the hydrate saturation increases, the connectivity between pores and throats shows a downward trend. Moreover, it is worth noting that pores and throats of model ST1 in the vertical and horizontal directions are basically not connected, when hydrate saturation is over 70%, and model ST08 and model ST06 have the same phenomenon, and the corresponding hydrate saturations are around 60% and 50%, respectively. This is because, as shown in Figure 3, the smaller the particle size, the smaller the radii of pores and throats will be, which makes pores and throats easier to be blocked. Figure 4 shows the hydrate saturations' effect on the average radius of pores and throats for three types of models under different particle sizes. And the intricate pore space was replaced by the spheres (pores) and cylindrical rods (throats). As can be seen in Figures 4(a) and 4(b), all the data could be utilized linear fitting and all the values of R 2 were approximately equal to 1, which means that the fitting was successful. It is apparent from Figure 4 that as the hydrate saturation increases, all the data decrease and the results of model ST1 are always the largest and that of model ST06 are always the smallest. This is because the particle size of ST1 is the largest and there are larger pores and throats; however, that of model ST06 is the complete opposite. Therefore, the larger the particle size, the larger the average pore radius and throat radius will be. Interestingly, it can be seen in Figures 4(a) and 4(b) that the smaller the particle size, the slower the decline rate will be. Specifically, the slopes of model ST1, model ST08, and model ST06 in Figure 4(a) are -3.56556, -2.77309, and -2.23594, respectively, while those in Figure 4(b) are -1.99065, -1.58598, and -1.26794, respectively. It appears that the larger the particle size, the more sensitive to the appearance of hydrate.
The fractal dimension is a crucial parameter to describe the pore structure's surface roughness, and the smoother the surface is, the smaller this parameter is. In this study, the box-counting method was utilized to calculate the fractal dimension, and the formula is as follows: where the N is the number of boxes which are covered with pores, r is the scale of the boxes, and D is the fractal dimension of the pore space. Figure 5 shows the particle sizes' effect on the fractal dimension for three types of models under different hydrate saturations. As shown in Figure 5, the fractal dimension of model ST06 is always the largest while that of ST1 is always the smallest. This is because the particle sizes from large to small are ST1, ST08, and ST06, and the smaller the particle size, the rougher the pore surface will be. It is important to note that the trends of these three models are highly similar. The fractal dimension of these three models all rises slightly to the maximum and then drops sharply as the hydrate saturations increase, and the value of the hydrate saturations corresponding to the highest point for model ST1, model ST08, and model ST06 is 30%, 27%, and 24%, respectively. This is primarily because when the hydrates start to generate, the initial size of them which are foreign matter to the surface of the pore is small. Therefore, at the beginning, the fractal dimension of three models will increase as hydrate saturations increase. However, when the hydrate saturation is high enough, the size of the hydrates will be about the same sizes as the particle sizes. At this time, the surface of the pores begins to become smooth, and the fractal dimension of the pores will decrease as the hydrate saturation increases.
The coordination number can reflect the connectivity of the pore system in the HBSs, and the smaller the coordination number, the worse the connectivity will be. Figure 6 shows the hydrate saturations' effect on the coordination number for three types of models under different particle sizes. It can be seen that as the hydrate saturation increases, the relative frequency of small coordination number gradually rises, and that of large coordination number gradually decreases. This is because of the presence of hydrates, some small throats are blocked, making the HBSs' connectivity deteriorate. It is worth noting that for the relative frequency of pores of which the coordination numbers are greater than 4 Geofluids 15, the smaller the particle size, the higher it will be. This is because the smaller particle size, the smaller throats will be, which will make the connectivity between pores and throats become better. Figure 7 shows the hydrate saturations' effect on the tortuosity for three types of models under different particle sizes. In Figure 7, there is a clear trend of increase in the tortuosity of HBSs with hydrate saturations increasing, and in the whole 5 Geofluids range of hydrate saturations, the tortuosity of model ST06 is always the largest while that of model ST1 is always the smallest. The reason for this phenomenon is that the small throats will be blocked due to the existence of the hydrates, making fluid flow through neighboring throats and resulting in the increase of the tortuosity. Moreover, there are more small throats in model ST06, but less small throats in model ST1. Therefore, the actual flow path of model ST06 is the most complicated, while that of model ST1 is the simplest which makes the tortuosity of model ST06 always the largest and that of model ST1 always the smallest. Interestingly, when the hydrate saturation is below 20%, the increasing rates of these three models are the fastest. This is because, at the beginning of hydrate formation, the number of the fine throats (200 μm-450 μm) in these models is large, which makes the throats more easily to be blocked by the hydrates.

Permeability Comparisons.
Permeability can influence the fluid flow capacity which will determine the gas production rate of the HBS [31]. This study utilized Darcy's law to calculate the HBS's absolute permeability through PNM and related parameters such as turbulence, compressibility, density, and the fluid's viscosity variability which were not considered. The formula is as follows: In the above formula, K h is the absolute permeability of the HBS, Q is the HBS's total flow, ΔP is the pressure gradient which is applied to the boundary, and in this study, this value was 2 MPa (input pressure was 12 MPa and output pressure was 10 MPa), μ is the fluid's viscosity and was set as 0.0011967 Pa·s which equals to water's viscosity at 11 MPa and 286 K, L is the fluid-flowing direction's length, and A is the fluid-flowing direction's cross-sectional area.
To better describe the relationship between the hydrate saturation, the particle size, and the permeability, this study utilized the Kozeny pore-filling model [32] and the Hou pore-filling model [33] to fit the related data. The specific formulas are as follows.
The Kozeny pore-filling model is as follows: where K 0 is the hydrate-free sediment's permeability, S h is the hydrate saturation, and a is the undetermined coefficient.
In model ST1, the value of a is -0.35954, in model ST08, this value is -0.27667, and in model ST06, this value is 0.14442. The Hou pore-filling model is as follows: where d is the undetermined coefficient, and in the model ST1, this value is 0.61178, in the model ST08, this value is 0.56952, and in the model ST06, this value is 0.82513. Figure 8 shows the hydrate saturations' effect on the permeability of the three types of HBSs under different particle sizes. It apparently shows that the Hou pore-filling model can fit data better than the Kozeny pore-filling model. This is because the larger the value of R 2 , the better fitting results will be, and the values of R 2 in model ST1, model ST08, and model ST06 in the Hou pore-filling model are larger than that in the Kozeny pore-filling model. As shown in Figure 8, with hydrate saturation increasing, the absolute permeabilities of three models begin to decrease fast and then decrease gradually. This is because hydrates will block the pores and throats which increase the resistance to flow. Moreover, at the beginning, the pores containing more gas and water make the formation rate of hydrates faster and the decline rate of permeability more obvious. In addition, as shown in Figure 8, the permeability of model ST1 is higher than that in model ST08, and that in model ST08 is higher than that in model ST06 which means the larger the particle size, the higher the permeability of the HBS will be. This is because the pores and throats in model ST1 are the largest, while those in model ST06 are the smallest, and the smaller the pores and throats, the greater the resistance for fluid to flow, and the lower the permeability will be [34].  Figure 5: The hydrate saturations' effect on the fractal dimension for three types of models under particle sizes.

Geofluids
3.3. Thermal Conductivity Comparisons. Temperature fluctuations have a great influence on the stability of HBS, so it is necessary to study the thermal conductivity of HBS. This part simulated the formation of hydrates under the gas saturation condition. And following Wu et al. [24], a partial differential equation based on Fourier's law considering multiple phases and multiscale homogenization theory was employed as follows:  8 Geofluids When the system reaches steady-state, the equation would be simplified to where the ðρc p Þ α is the heat capacity of the α phase, the ρ α is the density of the α phase, the c pα is the specific heat of the α phase, the λ α is the thermal conductivity of the α phase, and the T α is the temperature of the α phase. Moreover, in this part, the input temperature and output temperature were set as 300 K and 275 K, respectively, the thermal conductivities of sand, water, methane hydrate, and methane gas were set as 2.6 W/(m·K), 0.61 W/(m·K), 0.58 W/(m·K), and 0.03 W/(m·K), respectively. To better fit the related data, this part utilized the Dut model proposed by Wu et al. [24]: In the above formula, λ h is HBS's apparent thermal conductivity; t is the undetermined coefficient, and the value of t in the model ST1 is 3.38647, that in the model ST08 is 3.04341, that in the model ST06 is 2.60172; c is the undetermined coefficient, and the values of c in the model ST1, model ST08, and model ST06 are 0.00393, 0.00381, and 0.00378, respectively; λ 0 is the apparent thermal conductivity of hydrate-free sediment, and the values of λ 0 in model ST1, model ST08, and model ST06 are 0.69555, 0.73026, and 0.79023, respectively. Figure 9 shows the hydrate saturations' effect on the apparent thermal conductivity for three types of models under different particle sizes. It is apparent that as the hydrate saturation increases, the apparent thermal conductivities of these three models increase gradually. This is because the hydrate's thermal conductivity is higher than methane gas's thermal conductivity. Therefore, with increasing hydrate saturations, there are fewer and fewer gases in the HBS, making the HBS's thermal conductivity increase gradually. As shown in Figure 9, the apparent thermal conductivity of model ST06 is higher than that in model ST1 and model ST08, and the smaller the particle size, the higher the apparent thermal conductivity will be. This is because the smaller the particle size, the denser the HBS will be under the same condition, and the corresponding pores and throats will be smaller. Therefore, there are less gases in the pore space which make the HBS have a higher thermal conductivity.

Electrical Conductivity
Comparisons. Due to the great influence of hydrate saturation on the electrical conductivity of HBS, the electromagnetic remote-sensing technique is widely utilized to detect the distribution of hydrates. Therefore, it is important to understand the influencing factors of hydrate conductivity for developing a better exploitation plan. The apparent electrical conductivity is estimated using the following equation assuming no charge accumulation and no transient phenomenon following Wu et al. [24]: where the j ! is the current density, the j total is the total electrical flux going through the input face, the S is the area of the input face, the σ is the electrical conductivity of liquid phase, the V in and the V out are the input and output imposed potentials, and the L is the length of the specimen. And the formula for calculating the resistivity of HBS can be based on the previous study [24]. The specific formula expression is as follows: In the above formula, I is the HBS's resistivity index; R 0 is the hydrate-free sediment's electrical resistivity; R h is the HBS's electrical resistivity; e is the undetermined coefficient, and the values of e in the model ST1, model ST08, and model ST06 are 0.73732, 0.82165, and 0.92346, respectively; q is an undetermined coefficient, the value of q in model ST1 is 2.34337, that in model ST08 is 2.71112, and that in model ST06 is 3.51223; σ 0 is the hydrate-free sediment's apparent electrical conductivity; σ h is the HBS's apparent electrical conductivity. Figure 10 shows the hydrate saturation's effect on the electrical properties for three types of models under different particle sizes. It is apparent that as the hydrate saturations increase, the three types of models' resistivities rise rapidly, while the electrical conductivities show a completely opposite trend. This is because the electrical conductivity of liquid phase is higher than that of solid phase, and there will be less fluid in the pore space if hydrates in the pore space are more,  9 Geofluids making the electrical conductivities of the HBSs lower and the resistivities of the HBSs higher. In addition, the electrical conductivity of model ST1 is always higher than that of model ST08, and that of model ST08 is higher than that of model ST06, while the resistivities of HBSs are completely opposite. This is because the larger particle size is, the larger pores will be, and there is more liquid in the pore space which makes the electrical conductivity higher and the resistivity   Figure 9: The hydrate saturation's effect on the apparent thermal conductivity for three types of models under different particle sizes.
10 Geofluids lower. Besides, as shown in the image of electrical conductivity of HBSs, with hydrate saturations increasing, the electrical conductivity decline rate of model ST1 remains almost the same, while that of model ST08 and model ST06 gradually decreases, and the smaller particle sizes are, the obvious this phenomenon will be. This result suggests that the smaller particle sizes of the HBS, the much sensitive to electrical conductivity will be.

Conclusions
The physical characteristics of pore-filling HBS were systematically researched in this study. And the variation of pore characteristics' evolution, permeability, thermal conductivity, and electrical conductivity were studied under the different hydrate saturations and particle sizes. The conclusions are as follows: (1) The larger the particle sizes are, the larger the average radius of pores and throats will be, which makes the pores and throats will be more difficult to be blocked by the hydrates. For hydrate-bearing sediments, the connectivity of the pore system will get worse as hydrate saturations increase (2) The fractal dimension of the hydrate-bearing sediment will first increase and then decrease with increasing hydrate saturation. This is because only the size of the hydrates large enough will the pore surface become smooth. Moreover, the larger the particle sizes are, the smaller the fractal dimension will be (3) The smaller the particle sizes are, the larger the proportion of pores with large coordination number will be. This is because there are more small throats in the hydrate-bearing sediment with smaller particle sizes, which will make the connectivity between pores and throats become better (4) The tortuosity of hydrate-bearing sediment will increase gradually until to keep steady with increasing hydrate saturations. In addition, the larger the particle sizes are, the smaller the tortuosity will be and that of model ST1 is the smallest (5) With increasing hydrate saturations, the larger the particle sizes are, the higher the permeability of the hydrate-bearing sediment will be. This is because the smaller pores and throats, the resistance for fluid flowing will be greater. In addition, the values of R 2 in the Hou pore-filling model are greater than that in the Kozeny pore-filling model, which means the Hou pore-filling model fit the permeability data better (6) The apparent thermal conductivity of the hydratebearing sediment will increase with the hydrate saturation increasing. This is because the proportion of gas phase with lower thermal conductivity in pores gradually decreases. Moreover, the larger the particle sizes are, the higher the apparent thermal conductivity will be (7) As hydrate saturations increase, hydrate-bearing sediment's apparent electrical conductivity gradually decreases, while the resistivity index is opposite.

Geofluids
The smaller the particle sizes are, the lower the apparent electrical conductivity and the higher the resistivity index will be. Moreover, the larger the particle sizes are, the more linear the electrical conductivity decreases as hydrate saturations increase

Data Availability
The experimental data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.