Maximum Sizes of Fluid-Occupied Pores within Hydrate-Bearing Porous Media Composed of Different Host Particles

Hydraulic properties of hydrate-bearing sediments are largely affected by the maximum size of pores occupied by fluids. However, effects of host particle properties on the maximum size of fluid-occupied pores within hydrate-bearing sediments remain elusive, and differences in the maximum equivalent, incircle, and hydraulic diameters of fluid-occupied pores evolving with hydrate saturation have not been well understood. In this study, numerical simulations of grain-coating and pore-filling hydrate nucleation and growth within different artificial porous media are performed to quantify the maximum equivalent, incircle, and hydraulic diameters of fluid-occupied pores during hydrate formation, and how maximum diameters of fluid-occupied pores change with hydrate saturation is analyzed. Then, theoretical models of geometry factors for incircle and hydraulic diameters are proposed based on fractal theory, and variations of fluid-occupied pore shapes during hydrate formation are discussed. Results show that host particle properties have obvious effects on the intrinsic maximum diameters of fluid-occupied pores and introduce discrepancies in evolutions of the maximum pore diameters during hydrate formation. Pore-filling hydrates reduce the maximum incircle and hydraulic diameters of fluid-occupied pores much more significantly than grain-coating hydrates; however, hydrate pore habits have minor effects on the maximum equivalent diameter reduction. Shapes of fluid-occupied pores change little due to the presence of grain-coating hydrates, but pore-filling hydrates lead to much fibrous shapes of fluidoccupied pores.


Introduction
Natural gas hydrates vastly stored in marine sediments along the continental margins have a great potential to become one of global unconventional hydrocarbon energy resources [1][2][3]. Currently, the exploitation of this potential energy resource is still not economically feasible and requires innovative production methods and techniques [4,5]. New methods and techniques should be well evaluated by numerical simulators prior to field applications, and results of these numerical evaluations largely depend on proper characterizations of coupled processes and appropriate quantifications of physical properties of hydrate-bearing sediments [6][7][8][9].
Hydraulic properties (e.g., saturated water permeability, water retention curve, and relative permeability to water and gas) of hydrate-bearing sediments are inherently governed by pore-scale structures of the solid matrix and pores [17][18][19][20]. The pore space occupied by water and/or gas within hydrate-bearing sediments shrinks due to solid hydrate formation, and structures of fluid-occupied pores can be highly diverse due to different hydrate pore habits (e.g., graincoating and pore-filling) even though hydrate saturations are identical. These diverse pore structures are experimentally observed, and how they change during hydrate formation or dissociation has been quantified by using varieties of parameters with clear physical significances. Examples include porosity, shape factor, Euler characteristic of individual hydrate cluster, fractal dimensions, pore surface, and pore volume and size [21][22][23][24]. Various pore sizes (e.g., the critical, mean, and maximum pore sizes) have been correlated to the hydraulic permeability of porous media, and larger pore sizes generally lead to higher values of the hydraulic permeability [25][26][27]. Hydraulic properties of porous media are significantly affected by the maximum pore size [27][28][29], and the maximum pore size is experimentally and theoretically correlated with porosity, permeability, and particle size [30][31][32][33]. Grain sizes of marine sediments hosting natural gas hydrates are of a wide range from clays and silts to coarsegrained sands, and sand particle shapes are generally different [34][35][36]. Effects of host sediments properties (e.g., porosity, grain size, and shape) on the maximum size of fluidoccupied pores within hydrate-bearing sediments remain elusive, although papers have been published to clarify how the maximum size of fluid-occupied pores evolve with hydrate saturation during hydrate formation and dissociation [23,37].
Pores and fluid channels within most porous media in nature are generally nonspherical and noncircular [38][39][40]. For these irregularly shaped pores, several definitions of pore diameter are applied to quantify pore sizes. Examples include the equivalent diameter [41], the incircle diameter, and the hydraulic diameter [42]. The equivalent diameter λ e is a diameter of a circle having an area equal to the pore area. The incircle diameter λ i is determined by using the maximum ball method [43], and it is widely applied to pore network extractions from porous media [16,24,44]. The hydraulic diameter λ h is defined as λ h = 4A/P, where A is the cross-sectional area and P is the wetted perimeter of the cross-section. The hydraulic diameter is commonly used to simplify fluid flow in noncircular tubes and channels as round pipe flow [17,45,46]. For fluid-occupied pores within hydrate-bearing sediments, the maximum pore size is expected to decrease with increasing hydrate saturation. However, similarities and differences in the maximum equivalent, incircle, and hydraulic diameters of fluid-occupied pores evolving with hydrate saturation have not been well understood.
This study is aimed at clarifying the effects of host particle properties on the maximum size of fluid-occupied pores within hydrate-bearing porous media and further the understanding of different maximum pore diameters evolving with hydrate saturation. Grain-coating and pore-filling hydrates are randomly nucleated and grew within different artificial porous media to quantify the maximum equivalent, incircle, and hydraulic diameters of fluid-occupied pores at selected hydrate saturations, followed by analyses of maximum fluid-occupied pore diameters changing due to the presence of gas hydrates. Then, theoretical models for incircle and hydraulic diameter-related geometry factors are proposed based on the fractal theory, and these proposed models are further applied to provide insights into the hydrate saturation and morphology-dependent pore shape changes during hydrate formation.

Methods
Ten square images of artificial porous media shown in Figure 1(a) are generated by using the method of [47] for further numerical simulations of hydrate nucleation and growth in this study. Each of these porous media is constructed by randomly placing black particles with unrestricted overlap into a white square image until the desired porosity has been reached. The porosity is calculated as the ratio of white over total pixel numbers in a square image, and the square image has a side length of 300 pixels. The shape of a solid particle is characterized by using a concept of sphericity which is defined as S = 2 ffiffiffiffiffiffiffi ffi πA s p /P s , where A s is the area and P s is the perimeter of the solid particle [48]. Particle sphericity values are normally no bigger than 1, and S = 1 represents a circular solid particle. In addition, S = 0:86 stands for an elliptical solid particle with a major over minor diameter ratio of 2, and S = 0:76 for a diameter ratio of 3. Sphericity and size values of solid particles used for porous media constructions are summarized in Figure 1(b), and Figure 1(c) shows the intrinsic porosity and median particle diameter map for those artificial porous media.
Random nucleation and growth of grain-coating and pore-filling hydrates within artificial porous media is numerically simulated by using the method of [37]. These numerical simulations continually calculate minimum distances d s of fluid-occupied pore pixels to solid particle and hydrate pixels, and fluid-occupied pore pixels are selectively turned into hydrate pixels until the desired hydrate saturation S h has been reached. Different nucleation and growth preferences represent different hydrate pore habits. Grain-coating hydrate nucleation and growth is simulated by stochastically changing candidate pore pixels with d s = 1 px into hydrate pixels. Pore-filling hydrate nucleation and growth is modeled by preferentially seeding hydrate in candidate pore pixels with the highest d s , followed by randomly changing hydrate-touched pore pixels into hydrate pixels, and the value of the growth parameter is set to be 0.7 in this study. Each case of numerical simulations is 100 times performed to obtain probabilistically acceptable results. For more details, please refer to our previous paper [37].
The maximum equivalent Λ e , incircle Λ i , and hydraulic Λ h diameters of fluid-occupied pores within hydratebearing porous media are quantified from S h = 0 to S h = 0:8 at intervals of S h = 0:1. When a preselected hydrate saturation is reached, fluid-occupied pores are extracted from hydrate-bearing porous media by using the function named "bwlabel" with n = 4 (i.e., 4 connected pixels) in MATLAB 2016Ra, followed by calculations of the area and perimeter for all the fluid-occupied pores. Calculated values of the area and perimeter are further used to determine values of equivalent λ e and hydraulic λ h diameters, and the incircle diameter λ i is quantified by using minimum distance d s values. Assuming that there is a rectangular fluid-occupied pore with side lengths of 5 px and 6 px in hydrate-bearing porous media (Figure 2(a)), the equivalent diameter is calculated to be λ e = 6:2 px, the incircle diameter λ i = 5:0 px, and the hydraulic diameter λ h = 5:5 px according to their definitions ( Figure 2(b)). Values of the maximum equivalent, incircle, 2 Geofluids and hydraulic diameters of pores within hydrate-free porous media are determined and summarized in Figure 2(c) as intrinsic maximum pore diameters.

Results
The maximum equivalent Λ e , incircle Λ i , and hydraulic Λ h diameters of fluid-occupied pores within ten images of porous media containing gas hydrates are shown in Figure 3, and all the pore diameters decrease with increasing hydrate saturation S h . The average value of intrinsic maximum equivalent diameters Λ e0 of fluid-occupied pores within ten porous media (Figure 1(a)) is 154:7 ± 41:1 px with a confidence interval of 95.4% (the same below), and it decreases to 29:7 ± 16:7 px (19.2% of the intrinsic value) and 5:6 ± 1:3 px (3.6% of the intrinsic value) when hydrate saturation S h is 0.8 for grain-coating ( Figure 3(a)) and pore-filling ( Figure 3(b)) hydrates, respectively. The average value of intrinsic maximum incircle diameters Λ i0 of fluidoccupied pores within ten porous media is 47:8 ± 23:3 px, and it decreases to 28:0 ± 16:7 px (58.6% of the intrinsic value) for grain-coating hydrates ( Figure 3(c)) and 1:0 ± 0:0 px (2.1% of the intrinsic value) for pore-filling hydrates ( Figure 3(d)) when S h = 0:8. The average value of intrinsic maximum hydraulic diameters Λ h0 of fluid-occupied pores within ten porous media is 19:7 ± 9:1 px, and it decreases to 10:7 ± 3:7 px (54.3% of the intrinsic value) and 4:7 ± 0:6 px (23.9% of the intrinsic value) when S h = 0:8 for graincoating ( Figure 3(e)) and pore-filling ( Figure 3(f)) hydrates, respectively. It is obvious that pore-filling hydrates reduce values of the maximum equivalent, incircle, and hydraulic diameters more efficiently than grain-coating hydrates when S h = 0:8, and the maximum hydraulic diameter is the least sensitive to hydrate saturation. Host particle properties (e.g., intrinsic porosity, particle size, and sphericity) obviously affect the intrinsic maximum equivalent Λ e0 , incircle Λ i0 , and hydraulic Λ h0 diameters of fluid-occupied pores (Figure 2(c)) and introduce diversities (41.1 px for Λ e0 , 23.3 px for Λ i0 , and 9.1 px for Λ h0 ) in values of the intrinsic maximum pore diameter. All the intrinsic diversities decrease due to the presence of gas hydrates, and this implies that effects of host particle properties decrease 3 Geofluids with increasing hydrate saturation. In addition, pore-filling hydrates (Figures 3(b), 3(d), and 3(f)) have more significant effects on the intrinsic diversity reductions than grainfilling hydrates (Figures 3(a), 3(c), and 3(e)) when S h = 0:8.
In order to obtain further understandings of hydrate saturation and morphology-dependent maximum equivalent, incircle, and hydraulic diameters, Figure 4 shows normalized maximum equivalent Λ e * , incircle Λ i * , and hydraulic Λ h * diameters of fluid-occupied pores within hydrate-bearing porous media. Normalized maximum pore diameters are defined as ratios of maximum pore diameters within hydrate-bearing over hydrate-free porous media. All the values of normalized maximum equivalent diameters Λ e * of fluid-occupied pores within porous media containing grain-coating (Figure 4(a)) and pore-filling (Figure 4(b)) hydrates generally fall into the region between upper Λ e * = ffiffiffiffiffiffiffiffiffiffiffi ffi 1 − S h p and lower Λ e * = 1 − ffiffiffiffi ffi S h p curves, with Λ e * values partially lower than 1 − ffiffiffiffi ffi S h p when hydrate saturation is lower (i.e., S h < 0:5) for grain-coating hydrates and higher (i.e., S h > 0:5) for pore-filling hydrates. These two models are derived based on the assumption that gas hydrates uniformly grow into all pores with different sizes, and for their derivations, refer to [29,37]. Pore-filling hydrates reduce normalized maximum incircle (Figure 4(d)) and hydraulic (Figure 4(f)) diameters of fluid-occupied pores more significantly than grain-coating hydrates (Figures 4(c) and 4(e)). Normalized maximum incircle and hydraulic diameters of fluid-occupied pores decreasing due to the presence of grain-coating hydrates (Figures 4(c) and 4(e)) can be generally described by using Λ i,e * = ffiffiffiffiffiffiffiffiffiffiffi ffi 1 − S h p . Normalized maximum incircle diameter Λ i * of fluid-occupied pores within porous media containing pore-filling hydrates decreases with increasing hydrate saturation as Þ/2 when S h > 0 in a general trend (the blue dot curve in Figure 4(d)), and normalized maximum hydraulic diameter can be generally depicted by using [23,37] is applied to fit values of different normalized maximum pore diameters, and values of empirical parameters are summarized in Figure 4(g). The empirical model is an alternative form for the weighted average of theoretical models 1 − ffiffiffiffi ffi S h p for pore-filling hydrates and ffiffiffiffiffiffiffiffiffiffiffi ffi 1 − S h p for grain-coating hydrates [23]. Values of empirical parameters m and n are largely controlled by hydrate pore habits. For grain-coating hydrates, the normalized maximum equivalent diameter of fluid-occupied pores can be described by setting m = 0:067 and n = 70:6, normalized maximum incircle diameter by setting m = 0:71 and n = 5:98, and normalized maximum hydraulic diameter by setting m = 0:51 and n = 19:4. For pore-filling hydrates, the normalized maximum equivalent diameter of fluid-occupied pores can be described by setting m = 0:95 and n = 0:70, normalized maximum incircle diameter by setting m = 0:55 and n = 0:10, and normalized maximum hydraulic diameter by setting m = 0:50 and n = 0:79.

Discussion: Pore Shape Changes due to Hydrate Formation
Fractal theory is widely used to characterize pore-scale structures and investigate various physical (e.g., hydraulic and electrical) properties of porous media [49,50]. In these investigations, pores within porous media are treated as circles with equivalent areas in the two-dimensional space, and the maximum equivalent diameter Λ e can be calculated as [50] Λ e = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 π where A tot is the total area of porous media, ϕ is porosity, and D f is pore-size fractal dimension [28] which can be determined by using the box-counting method [51].  Figure 2: Equivalent, incircle, and hydraulic diameters of fluidoccupied pores within hydrate-bearing and hydrate-free porous media. (a) A pixel map of a rectangular pore occupied by fluids within hydrate-bearing porous media. Light gray pixels stand for the fluid-occupied pore with side lengths of 5 and 6 px, while dark gray pixels represent solid particle and hydrate pixels. Area A of the pore is 30 px 2 , and perimeter P is 22 px. Symbol d s represents the minimum distances of pore pixels to solid particle and hydrate pixels, and the maximum value ðd s Þ max in the rectangular pore is 3 px. (b) Values of equivalent λ e , incircle λ i , and hydraulic λ h diameters for the rectangular pore shown in (a). (c) Values of intrinsic maximum equivalent Λ e0 , incircle Λ i0 , and hydraulic Λ h0 diameters in hydrate-free artificial porous media shown in Figure 1(a).

Geofluids
For the case that size of a pore within porous media is quantified by using incircle diameter λ i (Figure 5(a)), area A p of the pore can be calculated as and c i is a geometry factor for incircle diameter. For references, the geometry factor for incircle diameter c i = 4/π for square-shaped pores, c i = ð3 ffiffi ffi 2 p Þ/π for regular triangleshaped pores, and c i = 1 for circular pores. Then, the maximum incircle diameter Λ i can be easily calculated as If hydraulic diameter λ h is used to quantify sizes of pores within porous media ( Figure 5(a)), area A p of an individual pore can be calculated as and c h is a geometry factor for hydraulic diameter. For references, the geometry factor for hydraulic diameter c h = 4/π for square-shaped pores, c h = ð3 ffiffi ffi 3 p Þ/π for regular triangleshaped pores, and c h = 1 for circular pores. Then, the maximum hydraulic diameter Λ h can be easily calculated as  Figure 3: The maximum equivalent Λ e , incircle Λ i , and hydraulic Λ h diameters of fluid-occupied pores within hydrate-bearing porous media evolving with hydrate saturation S h . The maximum equivalent diameter Λ e for grain-coating (a) and pore-filling (b) hydrates. The maximum incircle diameter Λ i for grain-coating (c) and pore-filling (d) hydrates. The maximum hydraulic diameter Λ h for grain-coating (e) and porefilling (f) hydrates.

Geofluids
If grain-coating hydrates uniformly grow in pores with different sizes and do not alter the shape of fluid-occupied pores ( Figure 5(b)), area A pc , incircle λ ic , and hydraulic λ hc diameters of these pores decrease with increasing hydrate saturation S h as A pc = A p0 ð1 − S h Þ, λ ic = λ i0 ffiffiffiffiffiffiffiffiffiffiffi ffi 1 − S h p , and λ hc = λ i0 ffiffiffiffiffiffiffiffiffiffiffi ffi 1 − S h p , respectively. In these equations, subscript 0   Figure 4: Normalized maximum equivalent Λ e * , incircle Λ i * , and hydraulic Λ h * diameters of fluid-occupied pores within hydrate-bearing porous media evolving with hydrate saturation S h . Normalized maximum equivalent diameter Λ e * for grain-coating (a) and pore-filling (b) hydrates. Normalized maximum incircle diameter Λ i * for grain-coating (c) and pore-filling (d) hydrates. The blue solid dot curve in (d) is drawn by Λ i * = ð1 − ffiffiffiffi ffi S h p Þ/2 when S h > 0 and Λ i * = 1 when S h = 0. Normalized maximum hydraulic diameter Λ h * for grain-coating (e) and pore-filling (f) hydrates. (g) A fitting model and its parameter values for descriptions of different pore diameters changing with hydrate saturation. The fitting model is drawn as black dashed curves in (a-f). 6 Geofluids represents the intrinsic (i.e., hydrate-free) condition and subscript c stands for grain-coating hydrates. Then, it is easy to obtain according to Equation (2), and based on Equation (4). Normalized geometry factor c ic * for incircle diameter in porous media containing grain-coating hydrates can be defined as according to Equation (6), and normalized geometry factor c hc * for hydraulic diameter as according to Equation (7). If pore-filling hydrates uniformly grow in pores with different sizes and the hydrate shape is identical with the intrinsic pore shape (Figure 5(c)), area A pf , incircle λ if , and hydraulic λ hf diameters of these pores decrease with increasing hydrate saturation S h as ffiffiffiffi ffi S h p Þ/2, and λ hc = λ i0 ð1 − ffiffiffiffi ffi S h p Þ, respectively. In these equations, subscript f represents pore-filling hydrates. Then, it is easy to obtain according to Equation (2), and based on Equation (4). Equation (10) is applicable for S h > 0, and c if = c i0 when S h = 0. Normalized geometry factor c if * for incircle diameter in porous media containing pore-filling hydrates can be defined as according to Equation (10), and normalized geometry factor c hf * for hydraulic diameter as according to Equation (11).
Values of the pore-size fractal dimension for fluidoccupied pores within porous media containing graincoating and pore-filling hydrates are summarized in Tables 1 and 2, respectively. The pore diameter ratio β defined as the ratio of the minimum over maximum pore diameters can be calculated by [52] and values of the pore diameter ratio within porous media containing grain-coating and pore-filling hydrates are shown in Figure 6. It is obvious that all the β values are generally smaller than 1:0 × 10 −2 , and the fractal theory can be used to analyze properties of porous media containing gas hydrates in this study [52]. Values of geometry factor c i for incircle diameter within hydrate-bearing porous media can be calculated by using Equations (3), where A tot = 300 × 300 px 2 . Geometry factors Gas hydrates uniformly grow in all pores with different sizes, and the shape of fluid-occupied pores before and after hydrate formation is unchanged. (c) Irregularly shaped pores in porous media containing pore-filling hydrates. Gas hydrates uniformly grow in all pores with different sizes, and the shape of solid hydrates is identical with that of intrinsic pores within porous media. Table  1: Pore-size fractal dimensions of fluid-occupied pores in porous media containing grain-coating hydrates.   The normalized geometry factor for incircle diameter changing due to the presence of grain-coating hydrates is shown in Figure 7(c) and pore-filling hydrates in Figure 7(d). It is shown that c ic * values generally stay close to the horizontal red line c ic * = 1 (Figure 7(c)), and c if * values go through an evolutionary process from below to above the red curve (Figure 7(d)). These discrepancies between numerical simulated c ic,f * data and corresponding theoretical models (i.e., Equations (8) and (12)) are mainly due to differences in hydrate pore habits since grain-coating and pore-filling hydrate growths can hardly follow the uniform and self-similar way ( Figure 5) strictly. Based on simulated data, an empirical model is proposed as  Figure 6: Pore diameter ratio β changing due to the presence of grain-coating (a) and pore-filling (b) hydrates within porous media. Regions colored in light green represent fractal theory applicable conditions (i.e., β < 1:0 × 10 −2 ).

Geofluids
shape changes in hydrate-bearing porous media during hydrate formation, while α = 2:8 and α = 5:1 set the lower and upper bounds. Geometry factors c h and normalized geometry factors c h * for hydraulic diameter evolving with hydrate saturation S h are shown in Figure 8. It is obvious that hydrate saturation and morphology effects on c i and c i * values are similar with those on c h and c h * values. Normalized geometry factor c hc * for hydraulic diameter in porous media containing grain-coating hydrates changes much more mildly compared with normalized geometry factor c hf * . The black curve in Figure 8(d) drawn by using Equation (13) agrees well with the general trend of c hf * values increasing due to the presence of pore-filling hydrates. Based on simulated data, another empirical model is proposed for the normalized geometry factor c if * prediction, which is with η = 15 describing the general trend, η = 3 setting the lower bound, and η = 50 setting the upper bound.
In general, grain-coating hydrates barely change while pore-filling hydrates significantly enhance c ic * and c hc * values, and perimeters of fluid-occupied pores decrease with an increase in grain-coating hydrate saturations but increase due to the presence of pore-filling hydrates according to their definitions. This implies that shapes of fluid-occupied pores do not obviously change due to the presence of grain-coating hydrates, but they alter to be more fiber-shaped when pore-filling hydrates occur.

Conclusions
This study numerically simulates random nucleation and growth of grain-coating and pore-filling hydrates within artificial two-dimensional porous media to quantify the maximum equivalent, incircle, and hydraulic diameters of fluidoccupied pores, and how maximum pore diameters evolve due to the presence of gas hydrates is analyzed. Theoretical and empirical models of geometry factors for incircle and hydraulic diameters of fluid-occupied pores are proposed based on fractal theory, and these proposed models together with simulated data are further applied to discuss effects of hydrate saturation and morphology on fluid-occupied pore shapes during hydrate formation. Main conclusions are drawn as follows. Intrinsic porosity, host particle size, and sphericity not only have obvious effects on the intrinsic maximum equivalent, incircle, and hydraulic diameters of fluid-occupied pores within hydrate-free porous media but also lead to discrepancies in the maximum pore diameters even through hydrate saturation and morphology are seemly identical.
Hydrate pore habits have relatively minor effects on the maximum equivalent diameter of fluid-occupied pores decreasing with increasing hydrate saturation. However, pore-filling hydrates reduce the maximum incircle and hydraulic diameters of fluid-occupied pores much more