Fundamental of swapping phenomena in naturally occurring gas hydrates

Amount of natural gas contained in the gas hydrate accumulations is twice that of all fossil fuel reserves currently available worldwide. The conventional oil and gas recovery technologies are not really suitable to gas hydrates because of their serious repercussions on geo-mechanical stability and seabed ecosystem. To address this challenge, the concept of methane-carbon dioxide (CH4-CO2) swapping has appeared. It has the potential in achieving safe and efficient recovery of natural gas, and sequestration of CO2. By this way, the energy generation from gas hydrates can become carbon neutral. This swapping phenomenon has not yet been elucidated at fundamental level. This work proposes a theoretical formulation to understand the physical insight into the transient swapping between natural gas and CO2 occurred under deep seabed and in permafrost. Addressing several practical concerns makes the model formulation novel and generalized enough in explaining the swapping phenomena at diverse geological conditions.

formation involving the exchange of multiple guest gases that are basically methane and CO 2 (pure or mixed) requires a generalized theory. Importantly, it should be versatile enough in predicting the transient state of the CH 4 -CO 2 and -mixed CO 2 exchange occurred within a porous cavity formed by pure/saline water molecules in the permeable sediments. There is no generalized model available to provide in-depth analysis of this multicomponent swapping process. With this research gap, this work proposes a theoretical formulation addressing various issues of practical importance for fundamental understanding of the naturally occurring gas hydrate phenomena in permafrost zones and under sea floor. Experimental data are used at diverse geological conditions to validate the formulation.

Results
Theory of gas swapping in hydrate lattice. A rigorous model is formulated for fundamental understanding of the transient sate of clathrate hydrates involved in the CH 4 -CO 2 (pure/mixed) swapping process. For this, the following chemical reaction is considered to represent the transition between gas, water and hydrates as H in which, g denotes the guest gas species, n H the hydration number, w the water species and h the hydrate component. Four phases may coexist there, namely gas (G), aqueous (A), solid ice (I) and hydrate (H). With this, the following practical aspects are taken into account to characterize the hydrate phenomena: • gas hydrates mostly occur in permeable marine sediments in presence of water (pure/with salt ions) 17 • hydrates are likely to form in the interstitial pore space between porous particles 18,19 • porous medium consists of irregular 3D particles with their uneven distribution • pores of these particles are further irregular in size and shape • these nanometer-sized pores also participate in hydrate formation and decay 1 . It is confirmed through the seismic survey studies conducted in the gas hydrate field of Alaska 20 , Blake Ridge 20,21 and Mackenzie Delta 22 • surface renewal is inevitable because of the barrier (i.e., solid hydrate film) grown during hydrate formation and decayed during dissociation at the interface 18 • pure carbon dioxide and methane, and 3-20 mol% CO 2 in air or with N 2 gas form sI hydrate structure that consists of two small 5 12 cages and six large 5 12 6 2 cages per unit cell 23 , which is evident through PXRD pattern 5 • pure and mixed CO 2 (N 2 major and CO 2 minor) effectively act as the replacement agent. In fact, mixed CO 2 has proven a better agent than pure CO 2 4 this is because CO 2 gets preference to occupy large cages, and in small cages, CH 4 and N 2 compete to each other for better occupancy 24 . Thus, our formulation attempts to accommodate both options.
• along with sand, glass beads are used as the second porous medium that mimics the effect of sediments for sand, sandstone and kaoline clays 25 .
Activity. Activity of water (a w ) has a great influence on hydrate dynamics. In the natural hydrate bearing atmosphere, it takes into account the collective effect of water in its pure form (subscript 'Pure'), and in presence of salt ion (subscript 'SI') and porous medium (subscript 'PM') as where, a w, Pure = γ w x w , in which, γ w represents the activity coefficient of water in water-gas mixture and x w the concentration of that water. The estimation of a w,Pure is briefly highlighted later. Further, the Pitzer model 26 is employed to formulate the water activity in electrolytic solution as where, MW is the molecular weight of water and m l the molality of solute species l (cation/anion/neutral). The osmotic coefficient, φ is estimated from 26 neutral ions in the system, respectively. The double summation indices, c < c′ and a < a′ refer to all the distinguishable pairs of dissimilar cations and anions, respectively. The molality of salt ions is obtained from: m l = n SI /W A , in which, n SI denotes the moles of salt ions and W A the amount of water remained in the aqueous phase, which is estimated from: W A = W in (1 − WC). Here, W in is the amount of water present prior to hydrate formation and WC the water conversion to hydrates during growth phase.
To formulate the activity of water in presence of porous medium, a couple of practical issues are to be addressed concerning the occurrence of gas hydrates in the irregular nanometer-sized pores with disordered capillaries of the distributed porous particles. Fractal theory 27 is used to model the activity, a w, PM in which, V w is the molar volume of water, k a constant, D f the fractal dimension of the pore edge, σ ∞ the interfacial energy between planar interfaces (0.0267 J.m −2 ) 28 , δ the Tolman length (0.4186 nm) 29 and R the universal gas constant. This modeling equation is applicable to both the irregular capillaries and pore, and thus the radius, r corresponds to both the hydrate core and pore. Here, we propose k as a linear function of r as: k = ar + b; a and b are the coefficients, values of which can be determined from the experimental data 29 .
Driving Force. The driving force is proposed in terms of chemical potential (μ) that considers the combined effect of temperature, pressure and composition 19 . It is expressed as where, μ w H and μ w L are the chemical potential of water in the filled hydrate and the liquid phase, respectively. Obviously, this Δμ is positive for hydrate formation and negative for dissociation.
Interstitial Reaction Surface Area. As stated, clathrate hydrates mostly form in the interstitial pore space between porous materials when small guest molecules (<0.9 nm) contact water at high pressure and low temperature. To formulate the cage dynamics, one needs to first find the total surface area of the irregular 3D particles 30 from: A = eV 2/3 , in which, V is the total volume of porous medium and e the scaling factor. For spherical shape, e is obviously 4.836 and for irregular solids, it is typically 6.2918. Now, the total volume of porous material is computed by subtracting the pore volume (V p ) from bulk volume (V b ) of the porous media with V b = m/ρ b , where m and ρ b denote the mass and bulk density of the porous material.
As forward reaction [in equation (1)] proceeds, the solid hydrate film starts increasing in size. Acting as a barrier at the interface, this film leads to decrease the contact area devoted for hydrate growth. Obviously, the reverse trend is true for backward reaction. Accordingly, the concept of effective surface area (A e ) is introduced as: A e = βA. As mentioned, the surface renewal is one of the important aspects of practical importance and it is proposed through updating the weighting factor as: β 0 = β 0 exp(Ct). Here, β 0 denotes the surface area adjustment factor, and the constant, C is negative for hydrate formation and positive for dissociation.
Reaction Kinetics. Prior to formulating the transient swapping process, one should kinetically model the hydrate formation and dissociation, and then couple them with certain flexibility in their movement. For this, apart from the driving force, Δμ, the hydrate kinetics is greatly influenced by the water consumption rate. With this, the reaction rate is formulated 12 as In which, Integrating and then simplifying, one obtains Here, α is a tuning parameter defined as the ratio of the highest amount of net guest gas consumed and the total amount of that gas ideally occupied in all cavities. This formulation will be used to represent the guest gas dynamics during hydrate formation and growth in salt water with porous media. Hydrate dissociation: The dissociation kinetics is governed by This is the modified form of equation (7), in which n H H O, 2 denotes the moles of water in hydrate phase. Substituting all relevant terms and simplifying, Here, n 0 denotes the total moles of guest gas present in the hydrate phase. This is the final formulation for hydrate dissociation in porous media and saline environment. Gas swapping: As stated, the swapping is a nondestructive process that proceeds with executing a dual mechanism of energy production and greenhouse gas sequestration. During this process, the pure/mixed CO 2 is introduced in the existing CH 4 hydrate bearing sediments at a reduced temperature and pressure to disturb the hydrate-gas equilibrium 32 . This leads to dissociating methane hydrates and forming mixed hydrates at the same time.
At this point, it should be noted that methane and nitrogen gas compete to occupy mostly small cages, while carbon dioxide preferentially occupies large cages without any challenge from other guests 24 . The thermodynamic mechanism for guest molecule exchange is typically controlled by the prevailing temperature, pressure and component composition 33 , on which, the chemical potential truly depends. The replacement continues through hydrate formation along with dissociation until the transition of equation (1) 14 . This thermodynamic control mechanism can also be explained by estimating the changes in free energy of gas hydrates 33 .
With this complicacy, we formulate the swapping dynamics in terms of methane displacement rate as The subscript f and d refer to the formation and dissociation of hydrate, respectively. Here, ξ is the ratio of fractional occupancy of CH 4 to that of replacement agent (subscript 'RA') (i.e., pure/mixed CO 2 ) in both small and large cages, which can be determined by the use of Langmuir type expression (equation (24)). Rewriting equation (14), Using integrating factor method, one can get the following form  This is the final theoretical representation of the proposed replacement kinetics that will be used in the sequel to analyze the transient swapping behavior of multicomponent gases with porous media and salt water.
Comparison to experiments. The replacement process in hydrate-bearing sediments is really complex because of the formation of CH 4 -CO 2 mixed hydrate and simultaneous dissociation of CH 4 hydrate in a system consisting of four phases (i.e., porous medium, hydrate, water and natural gas). To describe the complex physicochemical phenomena, we would like to use the proposed theory through validating it by the real-time data for two different sets of experimental arrangements (see supplementary information). All model parameters are reported in Table S1. With this, to quantify the performance of the developed formulation, we have used the absolute average relative deviation (AARD). CH 4 -Pure CO 2 Swapping with Aqueous Brine and Sand. First we study the CH 4 replacement by pure CO 2 (99.9 mol%) with the use of reported data 34 for experimental setup 1 (see supplementary information). Clay layers are practically impermeable owing to compact packing between constituent particles 4 , whereas sand layers, which contain a significant amount of gas hydrates, are considered to be practically permeable. Therefore, we consider sandy-sediment to simulate a real NGH environment. For this, the aqueous brine used has a salinity of 3.35 wt% (Na 2 SO 4 ) and the sediment is formed by 20-40 mesh quartz sands with a porosity of 38.7%.
The replacement reaction is performed in the reactor with 90.0 v% gas, 1.1 v% water and 8.9 v% hydrate saturation at stated operating pressure and temperature. During swapping, the natural gas gets replaced in hydrate cavities by the carbon dioxide gas. As a consequence, the CH 4 concentration increases in the gas phase and decreases in the hydrate phase. Naturally, the concentration dynamics of CO 2 is reverse in those two phases. This experiment leads to an about 45% of methane replacement by pure carbon dioxide. It is evident in Fig. 1 that the developed model shows a promising performance in predicting the swapping behavior between methane and carbon dioxide that is also confirmed through the AARD values.
Further, a direct comparison is made in Fig. 1 between CH 4 swapping kinetics with CO 2 (simulation and experiment) and CO 2 /N 2 mixture (simulation only, since experimental data are not available). It is a fact that the density of injected CO 2 is high, while its permeability is reasonably low. Adding N 2 to CO 2 leads to increase the gas permeability. Moreover, blocking of flow channels owing to the formation of new hydrate from injected CO 2 /N 2 is less compared to injected pure CO 2 35 . Mainly because of these reasons, the CO 2 composition in hydrate phase is higher in case of CO 2 /N 2 mixed guest gas compared to pure CO 2 , which is quite obvious in Fig. 1 with respect to time.

CH 4 -Mixed CO 2 Swapping with Glass Beads.
In naturally occurring gas hydrate sites, injected carbon dioxide gas may transform to a liquid state at harsh conditions 5 . This leads to make injection and diffusion very unstable, yielding a low replacement rate and recovery. To overcome such weakness, the use of CO 2 /N 2 gas mixture is proposed and successfully performed in a field production test on the Alaska North Slope in 2012 5 . This mixed gas improves the replacement efficiency compared to pure CO 2 mainly because of replacing natural gas by CO 2 in large cages and that by N 2 in small cages 24 , along with the reasons mentioned in the last case study. Moreover, the formation condition of CO 2 hydrate is known to be more stable than that of CH 4 hydrate 31 , indicating the suitability of swapping process for long term CO 2 storage 24 . Keeping these issues in mind, the proposed theory is tested here with the mixed CO 2 gas as well in presence of glass beads that mimic the effect of sediments for sand, sandstone and kaoline clays, which have little to no swelling in contact with water 25 .
With this, the performance of the developed formulation is evaluated with simulating the compositional dynamics at a distance of 0.7 m from the inlet of the 8 m long reactor. Here, the CH 4 present in hydrate cavities is attempted to swap by injecting the replacement gas mixture (i.e., 20 mol% CO 2 and 80 mol% N 2 ) into the reactor at a rate of 100 sccm (Fig. 2a) and 200 sccm (Fig. 2b). The replacement reaction starts in the reactor with 73.67 v% gas, 14.6 v% water and 11.73 v% hydrate saturation for 100 sccm case, and 73.64 v% gas, 14.6 v% water and 11.76 v% hydrate saturation for 200 sccm case. The reactor contains a porous bed formed with glass beads.
The guest gas (i.e., CO 2 /N 2 mixture) injection rate truly affects the replacement kinetics and total run time 4 . For example, increasing gas injection rate leads to shallow penetration into the gas hydrates layers. It is because of short residence time with high flux, and thereby low replacement efficiency. On the other hand, the low flux of the guest gas leads to deep penetration and more extended stay, which results in improved replacement efficiency. Here, we have demonstrated the model prediction for two different gas mixture injection rates (i.e., 100 and 200 sccm). Due to the stated reason and as experimentally observed 4 , the 200 sccm case (Fig. 2b) leads to decrease the replacement efficiency compared to the case with 100 sccm injection rate (Fig. 2a). Further, the run time decreases with increasing injection rate as the methane gets depleted rapidly in the reactor 4 . It is obvious that the proposed model provides a reasonably good prediction of swapping data.
Injecting guest gas (e.g., CO 2 ) into in-situ natural gas hydrates in sediments leads to conversion over to CO 2 and mixed CO 2 /CH 4 hydrates. Now, this conversion is governed by the two primary mechanisms 35 . They are based on: (i) direct solid state conversion, and (ii) formation of new hydrates from injected guest gas and free water. In the first mechanism, the CH 4 hydrate directly converts into pure CO 2 and mixed CO 2 /CH 4 hydrate. This is a very slow process because of slow mass transport through hydrate and it dominates only when there is no sufficient free water available. The previous case study (Fig. 1) deals with only 1.1 v% water saturation, and thus the first mechanism dominates there, which leads to reach about 20% CO 2 composition in hydrate phase taking a reasonably long time (around 60 hr).
As far as the second mechanism is concerned, the injected CO 2 gas reacts with free water in the porous media and forms new CO 2 hydrate. This reaction is exothermic in nature and thus, it releases heat, which in turn further dissociates the surrounding CH 4 hydrate. Then, with the generated water, the CO 2 forms more hydrates. Note that this mechanism dominates if sufficient amount of free water is available there and it is faster than the direct solid state exchange mechanism. The CO 2 gas consumption from the mixed CO 2 /N 2 guest gas in Fig. 2 follows this mechanism (14.6 v% water saturation) and thus, it reaches about 20% CO 2 composition in around 10 hr, which is about 60 hr in the previous case (1.1 v% water saturation) (Fig. 1).
In the subsequent study, the guest exchange ability of CO 2 /N 2 mixture is shown at the prevailing condition by attacking and replacing the encaged CH 4 in hydrate bearing sediments. The comparative gas composition profiles are produced as a function of length for an injection rate of 100 sccm, keeping all other conditions same with the last test (i.e., pressure = 9.8 MPa, temperature = 275.15 K, replacement agent having 20 mol% CO 2 and 80 mol% N 2 , 100 μm average size glass beads, and 14.6 v% water saturation). With this, Fig. 3 compares the proposed model predictions with reference to the experimental data 4 for three different sampling ports situated at (a) 2.4 m, (b) 4 m and (c) 5.6 m in the 8 m long one-dimensional tubular reactor. Expectedly, the CH 4 composition decreases with time, and for CO 2 and N 2 , the trend is just reverse. As the guest gas is injected, the CH 4 is displaced, triggering spontaneous swapping reaction in the hydrate phase.
This test further proves that the proposed formulation is capable of predicting the replacement phenomena with reasonable accuracy that is indicated through the AARD values. It should be noted that for establishing the flow-front of gases in the sandy gas hydrate layer and a reliable estimation of natural gas recovery, this integrated dynamic compositional information along the length is inevitable 4 .

Discussion
Gas swapping inside the hydrate lattice is proposed as a nondestructive process that can provide the next generation energy for a couple of centuries and reduce the emission of greenhouse gas arising from anthropogenic activities. Addressing several fundamental features of gas hydrates and their occurrence in permafrost and under deep seabed at diverse geological conditions in theory, attempt is made to make the model realistic and generalized enough. Two porous media composed of unevenly distributed sand particles and glass beads are separately used to investigate the model performance in presence of pure water and aqueous brine. There is a key finding reported through the experimental investigation that CH 4 -CO 2 /N 2 provides a higher replacement efficiency than CH 4 -CO 2 4 . Therefore, the proposed formulation is tested with binary and multicomponent guest gases to show its versatility. With this, the developed formulation describes the swapping phenomena precisely and provides a promising performance with a close agreement with the experimental data.
This formulation can further be used • to understand the hydrate characteristics involved in natural gas storage and transportation, desalination, gas separation, and CO 2 capture and sequestration, among others • as a scale-up model for all the above-mentioned hydrate specific processes Providing clearer insight into how hydrates are generated, deposited and decayed, this model can play a crucial role in the assessment of total gas reserves in hydrate deposits that is still a highly volatile issue, involving a wide range of uncertainty.

Methods
Activity. Finding a w,Pure . First of all, we need to determine x w from: x w = 1−x g , when γ w for pure water is unity. Then one can obtain the composition of guest gas (x g ) by knowing the molality of that gas (m g ) from For the guest gas, μ L(0) and μ V(0) are the standard chemical potential in the liquid and vapor phase, respectively, and φ the fugacity coefficient that is obtained with the use of Soave-Redlich-Kwong (SRK) equation of state 36 . Here, μ V(0) is adopted as zero 37 , and μ L(0) is considered as a function of operating temperature and pressure 37 .
Finding a w,SI . The parameters, namely A φ , C φ and ψ, are function of temperature (T) as whereas, λ and ζ are dependent on temperature (T) and pressure (P) as Values of constants, c 1 to c 6 and C 1 to C 6 , for all these five parameters are available in literature 26 . Further, the ionic strength (I) is defined as: = ∑ I mz l l l 1 2 2 , where z l is the charge of ion l. Moreover, equation (4) includes Z that is expressed as: = ∑ Z m z l l l . Second virial coefficient, φ B ca has the following representation where, β (0) , β (1) and β (2) are the temperature dependent parameters and equation (18) is used to calculate them. Note that α ca is equal to 2.0 for univalent and 1.4 for higher valence pairs 26 . Besides, an another second virial coefficient, Φ that accounts for the interaction between the ions of equal sign (l−m), is given as Here, Θ is a single parameter for each pair of anions or cations, and the functions, E Θ lm (I) and E Θ′ lm (I), take into account the electrostatic unsymmetrical mixing effect and depend on ionic strength and the type of electrolyte pair. The single electrolyte third virial coefficient, C ca is estimated as: Finding a w,PM . It is governed by the following equation where, ΔP is the difference in pressure between the liquid and hydrate phase. The surface effects of the pore edge of irregular capillaries and irregular pores on ΔP are taken care of through: , in which, p is the perimeter (=2πkr Df ), s the area (=πr 2 ) and θ the wetting angle between hydrate and pore wall, which is zero 29 . where, μ w 0 is the chemical potential of water in empty hydrate cages, μ Δ w H the difference in chemical potential of water between empty and filled hydrate cages in hydrate phase, which is calculated from 31 . Like μ w H , μ w L includes μ w 0 and μ Δ w L according to equation (23). This μ Δ w L , which represents the difference in chemical potential of water between the empty hydrate cages in hydrate phase and the liquid phase, is estimated from the standard chemical potential difference of water for gas hydrate at reference temperature and absolute zero pressure, Δ V w L the difference between molar volume of water in hydrate and liquid phase, and Δ h w L and Δ C p L the enthalpy and heat capacity difference between empty hydrate cages and liquid water, respectively.