Formulating formation mechanism of natural gas hydrates

A large amount of energy, perhaps twice the total amount of all other hydrocarbon reserves combined, is trapped within gas hydrate deposits. Despite emerging as a potential energy source for the world over the next several hundred years and one of the key factors in causing future climate change, gas hydrate is poorly known in terms of its formation mechanism. To address this issue, a mathematical formulation is proposed in the form of a model to represent the physical insight into the process of hydrate growth that occurs on the surface and in the irregular nanometer-sized pores of the distributed porous particles. To evaluate the versatility of this rigorous model, the experimental data is used for methane (CH4) and carbon dioxide (CO2) hydrates grown in different porous media with a wide range of considerations.

and it predominates in the Earth's natural environments; and cubic structure II usually forms with larger (0.6-0.7 nm) guests in the man-made environments 10 .
Despite worldwide abundance of gas hydrate and its potential as an energy source of the future, the growth kinetics of this crystalline structure remains poorly understood 19 . To address this issue, in this study, a mathematical formulation is developed in the form of a physical model that is capable and versatile enough in precisely explaining the formation kinetics of various gas hydrates. Perhaps, there is no such generalized formulation reported in literature to predict the real-time growth behavior.

Results
Formulating hydrate formation kinetics. Hydrate formation seems to occur at the interface between the bulk guest and aqueous phases. This formation may be driven by the difference existed in temperature 20,21 , pressure 22 , gas composition 23 or in fugacity [24][25][26][27][28] . Here, the driving force is proposed as the difference in chemical potential of water in the aqueous phase ( µ ∆ w A ) and water in the hydrate phase ( µ ∆ w H ). Clathrate hydrates mostly form in the interstitial pore space between porous particles 29 . Accordingly, the consumption of guest gas is proposed to vary proportionally with the said chemical potential difference, along with the water transformation rate and the total particle surface area (A). This yields, the residual mole of water, K 0 the intrinsic rate constant, ΔE a the activation energy, T the temperature and R the universal gas constant (8.314 J. mol −1 .K −1 ). This modeling equation is formulated by assuming the first-order reaction kinetics for water transformation in terms of n H O,A 2 . Further, it considers that the kinetic constant (K) is temperature dependent and represented by an Arrhenius-type equation as: (2) 0 a As time progresses, the hydrate film increases in size and it acts as a barrier at the interface. This leads to decrease the contact area involved in hydrate growth between the bulk guest and aqueous phase. Accordingly, the concept of effective surface area (A e ) is introduced 29 as: The surface area adjustment factor, β lies between 0 and 1. Actually, β needs to be tuned by the use of an optimization technique. Now, one needs to replace n H O,A 2 by the mole of guest gas (n gg ), for which, the following equation can be used: This is obtained by using the following correlations:  Here, t denotes the time and α an adjustable parameter that is defined as the ratio of highest value of the net amount of guest gas consumed during the hydrate formation and the total amount of that gas ideally occupied in all cavities.
Targeting to the formation of CH 4 and CO 2 gas hydrates, both of which are sI type hydrate, one can compute µ w H from 30 : the standard chemical potential difference of water for gas hydrate at reference temperature and absolute zero pressure, and it is adopted as 1202 J.mol −1 (ref. 31), ∆V w A the difference between molar volume of water in hydrate and aqueous phase, ∆h w A the enthalpy difference between empty hydrate lattice and liquid water, and a w the activity of water.
Along with the interstitial pore space between porous materials, gas hydrates are likely to form and grow inside the nanometer-sized pores of those materials. In this regard, one can see the experimental evidence provided in ref. 32. To formulate this growth kinetics, it is quite realistic to consider that the pores are irregular in shape and size (schematic in Fig. 1), and they are present in the distributed particles. Further, it is considered that a thin bound water monolayer having a thickness of 0.4 nm 33 is present on the pore wall. With these, assuming self-similar characteristics of the pore edge existed in the hydrate media, the following form of expression is used to compute the activity of water (a w ) for both the growth sites as: where, γ w represents the activity coefficient of water (assumed unity) 31 , x w the composition of water, V w the molar volume of water, D f the fractal dimension of the pore edge, r core the radius of hydrate core, σ ∞ the interfacial energy between planar interfaces and δ the Tolman length. Note that k is a linear function of pore radius (r pore ) as shown later. Now, one can simplify the case with considering hydrate formation in the regular pores of the porous materials and on their effective surface. Supposing cylindrical pores and circular pore edge, the following expression 30, 34 is used to estimate the activity of water (a w ) as: For this, the available experimental data sets are used and the model performance is quantified in terms of the absolute average relative deviation (AARD) that is expressed as: where, n dp is the total number of experimental data points, and WC e and WC p are the experimental and model predicted water conversion to hydrates (%), respectively. This WC is estimated as: The subscript gg, H 2 O and H refer to the guest gas, water and hydrate, respectively. It should be noted that the conditions used in modeling and experiments in references are same.
Silica gel. With silica gel and a binary gas mixture of CO 2 -N 2 (17-83%), the following four versions of the model are proposed to find the best performing case: Case I: Only CO 2 forms hydrate on the surface and in the irregular pores of the distributed porous particles. Case II: Both CO 2 and N 2 form hydrate on the surface and in the irregular pores of the distributed porous particles. Case III: Only CO 2 forms hydrate on the surface and in the regular pores of the distributed porous particles. Case IV: Both CO 2 and N 2 form hydrate on the surface and in the regular pores of the distributed porous particles. Figure 2a and b compare the above four versions of the kinetic model with reference to the experimental data 35 in terms of water conversion to gas hydrates formed in silica gel at two different operating pressures. In this study, the distributed porous particles are used at their arithmetic mean radius. The model parameters, namely α, β, ΔE and K 0 are optimized using the generalized reduced gradient (GRG) method, and they are reported in Table 1. With this, the best performance is achieved by Case I of the developed model as evident from the AARD value (7.94% in Fig. 2a and 3.91% in Fig. 2b) mainly because of considering irregular pores for hydrate formation and growth. Further, comparing two cases (I-II or III-IV) it becomes obvious that N 2 does not have any significant role in hydrate formation. In the sequel, thus the best performing version (i.e., Case I) of this model will further be investigated to gain insight into the formation kinetics.
In the next study, attempt is made to investigate the effect of particle size on water conversion to the hydrate phase. Determining optimal parameter values ( Table 2) and then using them, it is investigated in Fig. 3 that the developed kinetic model performs closely, despite a difference existed in particle size. This is achieved by selecting an optimal value for β and ΔE a against each particle size as shown in Table 2. This scope of retuning leads to make the model almost unaffected by the particle size as indicated in terms of AARD values highlighted in Fig. 3 (~7.93% in all three cases).
Further, Fig. 4 evaluates the formulation made in the form of a model for hydrate formation using the experimental data taken from literature 36,37 . In these experiments, a variation is made in the feed gas along with the operating pressure. For the binary feed consisting of CO 2 and H 2 (Fig. 4a) and the ternary feed of CO 2 , H 2 and C 3 H 8 (Fig. 4b and c), only CO 2 forms gas hydrate as indicated before. Along with following the real-time growth trend, the model shows its promising performance in the aspect of the degree of closeness achieved between the  (Fig. 2b). Moreover, the runs are performed with a feed gas mixture of CO 2 -N 2 (17-83%) in the bed of silica gel with 45-75 µm particle size distribution having a pore radius of 30 nm. The percent AARD values are given in the figure against each of the four cases.
predicted and experimental figures. This is reflected through the error calculated in terms of AARD (i.e., 7.94% in Fig. 4a, 4.99% in Fig. 4b and 4.94% in Fig. 4c).
Silica sand. Now attempt is made to test the model for CO 2 hydrate formation in a different porous media, namely silica sand. Unlike the silica gel, it has reasonably small pores. Here, two sets of results are produced in Fig. 5a and b, which are different in terms of their operating pressures. Finding the optimal parameter sets ( Table 1) and using them, it is evident that the proposed formulation is capable enough in predicting the real-time formation behavior. In this regard, one can see the AARD values provided in the figure itself.
It is observed from the experimental investigation 37 that to a certain extent, the initial rate of water conversion increases with the increase of pressure. The subsequent parts of the conversion rate profiles follow the similar trend but with a larger magnitude at a higher pressure. Now, one can closely observe that the real-time water conversion at 4.5 MPa (Fig. 5a) is relatively slow at the beginning and then it picks up speed. This typical behavior (sigmoidal shape) can be approximated by the response of a first-order system with a time lag. On the other      (Fig. 5a) and 5.5 MPa (Fig. 5b).
hand, the water conversion at 5.5 MPa (Fig. 5b) leads to the response of a first-order system only (negligible time lag). Now, the governing equation in the proposed model [i.e., Equation (7)] is quite similar to the model of a first-order system (with no time lag) and thus, the model predicts the conversion rate at 5.5 MPa better than that at 4.5 MPa. However, this difference in initial conversion rate between 4.5 and 5.5 MPa does not exist for silica gel as shown in Fig. 4b and c before, and thus, the AARD values are quite close between them. It is fairly true that the natural hydrate formation process continues for a long period of time. Keeping this issue in mind, the model predictability is further tested in Fig. 6 for continuously about 120 hour (5 days). In case of pure CO 2 gas (99.9%), the experimental data are taken from literature 38 . It is evident that the proposed model shows a good prediction with a reasonably low AARD of 6.82%.

CH 4 Hydrate.
To evaluate the developed formulation and its versatility, further the formation and growth of CH 4 hydrate are considered in presence of two porous media, namely silica sand and hollow silica. As stated earlier, the same experimental conditions are used in the model simulation.
Silica sand. Figure 7a depicts the performance of the kinetic model with reference to the experimental data 1 in the aspect of water conversion to CH 4 hydrate in silica sand. The pure CH 4 gas is used for the hydrate formation and subsequent growth. This study is performed at 8 MPa and 277.15 K with the average particle size chosen for the distributed range of 560-1300 µm. The pore size is considered as 0.9 nm. Using the identified model parameters (Table 1), the model shows an excellent agreement with the data with an AARD of about 4%. This is achieved by addressing a couple of practical issues in the model formulation as stated before. It should be noted that this test is conducted for a long period of time (about 4200 min ( = 70 hour)) typically involved in the natural hydrate formation process.
Hollow silica. To predict the real-time kinetic profile of CH 4 hydrate growth, the formulation made in this study is used for the hollow silica distributed in the range of 30-70 µm. The system operates at 8 MPa and 278.2 K. Taking average particle size, the ratio of the mass of hollow silica and the volume of water is considered as 1:4, 1:6 and 1:8 in Fig. 7b,c and d, respectively. From the results, one can see the promising performance of the model, indicating that the formulation has well taken the issues that have practical relevance.

Discussion
Here, a physical model is developed to understand the hydrate formation phenomena. This model is novel in that it considers the clathrate hydrate formation in both the interstitial pore space between porous materials and inside the nanometer-sized pores of those materials. By considering chemical potential as a driving force for growth, the combined effect of temperature, pressure and composition is taken into account. More importantly, the proposed formulation addresses a couple of practical issues concerning irregularity in the pores of distributed particles, surface tension effect in the pores, among others. Excellent agreement is achieved between the model prediction and experimental data for several porous media, and this is also reflected through the AARD values. Thus, it can be concluded that the proposed formulation is rigorous and versatile enough to represent a generalized model in predicting the formation kinetics of clathrate hydrates. This model can further be improved by renewing the surface area of the water in contact with the hydrate gas with time, and considering pore size distribution in the non-spherical porous particles.

Methods
Estimating µ w A . To use Equation (9) that models µ w A , one needs to estimate ∆h w A , for which, the following form is recommended, Figure 6. Performance evaluation with CO 2 hydrate formation in silica sand. For this study, the data sets 38 are available at 277.2 K in presence of the said porous medium having particles distributed in the range of 100-500 µm with a pore size of 0.9 nm. A pure CO 2 (99.9%) is used in the experiment conducted at 3.5 MPa.
To model the activity of water (a w ) in a porous medium, the following form 30 is used: where, ΔP denotes the difference in pressure between the aqueous and hydrate phase. x w is determined based on the guest gas present in aqueous solution (x gg ) as: Here, µ gg A and µ gg V are the chemical potential of guest gas in the aqueous and vapour phase, respectively. Now one can calculate x gg by estimating the molality of the respective guest gas (m gg ) from the following equation: where, m w is the number of moles of water per kg of water (55.56 mol/kg) 31 , y gg the mole fraction of guest gas in vapour phase, µ A gg (0) and µ V gg (0) the standard chemical potential of a guest gas in the aqueous and vapour phase, respectively, φ gg the fugacity coefficient of guest gas that is estimated here by using the Soave-Redlich-Kwong equation of state 39 , and γ gg the activity coefficient of guest gas (assumed as one) 40 . Note that the µ V gg (0) is adopted as zero 31 .
The µ A gg (0) is a function of operating temperature and pressure, and it is calculated by the following equation for methane hydrate 31 : in which, C 1 to C 10 are the coefficients. Similar form of equations is also reported for other gas hydrates 40,41 .
As far as ΔP is concerned, it is proposed to consider the growth of hydrate in irregular nanometer-sized pores. For this, the following form 33 is used: H A in which, σ H−A is the surface tension of water between aqueous and hydrate phase, and θ the contact angle between water and porous media, which is zero. Here, L and S are the perimeter and area of the pore edge for the irregular capillaries, respectively. According to the fractal theory 42 , L can be calculated using the fractal dimension of the pore edge (D f ) and pore radius (r pore ) as: where, k is a linear function of r pore , pore in which, a and b are the coefficients that can be estimated from the experimental data 33 . Here, it is supposed 33 that the area of hydrate core, S is same as that of circular hydrate core. Similarly, the perimeter of hydrate core (l) (Fig. 1b) can be expressed in terms of r core as: Further simplifying and rearranging, ΔP can be expressed as: in which, σ H−A is represented as follows 43 : where, σ ∞ is adopted as 0.0267 J.m −2 (ref. 44). The thickness of an interfacial region between solid (ice) and aqueous (water) phase (δ) is commonly referred as Tolman length, which is equal to 0.4186 nm 45 . In addition, the solid-liquid interfacial curvature (κ) is considered as a function of r core and D f , and it is given as 33 : core 2 D f Substituting Equations (25)(26)(27) in (15), one obtains the model of a w represented in Equation (10).
To simplify the case with hydrate formation in the regular pores and on the effective surface of the porous particles, it assumes cylindrical pores and circular pore edge, for which 30 in which, C ij represents the Langmuir constant of gas component j in an i type cavity. The fugacity of gas component j in the hydrate phase (f j ) is estimated from the Soave-Redlich-Kwong (SRK) equation of state. This f j is assumed same with the fugacity of component j in the gas phase 46 . Now, C ij is computed from 30  Here, K denotes the Boltzmann's constant, R the cell radius of hydrate and ω(r) the spherically symmetric cell potential, which is obtained from 47 : a r a ( ) 4 2 2 2 2 (31) 12 6 where, the constants 47 , ε, σ and a denote the maximum attractive potential, the cores distance at zero potential and the radius of the spherical core, respectively.
Estimating particle surface area for hydrate growth. Along with the pores, as stated, the surface of the porous particle is also involved in hydrate formation and growth. In this light, the concept of effective surface area (A e ) is introduced in the proposed model. Now the total surface area of the porous material (A) is estimated by multiplying the individual particle surface area (A pi ) with the total number of particles (n tp ) present in the bed as: tp pi Usually, the size of the porous material is known in terms of its diameter (d p ) and hence, it is easy to find the A pi assuming spherical particle from: On the other hand, n tp is obtained from: where, V tp denotes the total volume of the porous media and V pi the volume of a single particle. The V tp is calculated by subtracting the volume of water required to fully saturate the fixed bed (V ws ) from the total volume of the bed (V b ). Knowing V ws , one needs to determine V b from: Here, d b and h b are the diameter and height of the fixed bed of a porous medium.
Data availability. The data sets that support the findings of this work are available from the corresponding author upon request.