Modeling recovery of natural gas from hydrate reservoirs with carbon dioxide sequestration: Validation with Iġnik Sikumi field data

Fundamental understanding of guest gas replacement in hydrate reservoirs is crucial for the enhanced recovery of natural gas and carbon dioxide (CO2) sequestration. To gain physical insight into this exchange process, this work aims at developing and validating a clathrate hydrate model for gas replacement. Most of the practical concerns associated with naturally occurring gas hydrates, including hydrate formation and dissociation in interstitial pore space between distributed sand particles in the presence of salt ions and in irregular nanometer-sized pores of those particles, irregularity in size of particles and shape of their pores, interphase dynamics during hydrate formation and decay, and effect of surface tension, are addressed. An online parameter identification technique is devised for automatic tuning of model parameters in the field. This model is employed to predict the laboratory-scale data for methane hydrate formation and decomposition. Subsequently, the model is validated with the field data of the Prudhoe Bay Unit on the Alaska North Slope during 2011 and 2012. In this Iġnik Sikumi field experiment, mixed CO2 (i.e., CO2 + N2) is used as a replacement agent for natural gas recovery. It is observed that the proposed formulation secures a promising performance with a maximum absolute average relative deviation (AARD) of about 2.83% for CH4, which is even lower, 0.84% for CO2 and 1.67% for N2.


Motivation and objective
Naturally occurring hydrates form once an appropriate guest gas like methane comes in contact with the water molecules at a reasonably high pressure (i.e., more than the atmospheric pressure), and/or the low temperature 26 . The most common guest gases such as methane and carbon dioxide are recognized as s-I hydrate formers 5 , and their unit cell comprises of two small and six large cavities made from the hydrogen-bonded water molecules network, as depicted in Fig. 1. The natural gas hydrates (NGH) are truly a great source of hydrocarbon energy. The potential energy content is about twice that of all the other fossil fuels combined together 1 . Deep sea-beds and permanently frozen grounds are the favourite and habitable hosting sites for the NGH. In these reservoirs, usually the hydrates are enclosed in the unconsolidated porous sand layers/sediments, and accompanied by the saline water. Among these, former element (porous medium) acts in favour of the hydrate formation, and later one (seawater) typically hinders their growth. Thus, it is inevitable to overlook their influence on the naturally occurring gas hydrates.
The gas hydrates prefer to cultivate themselves in between the interstitial pore spaces as well as in the internal nano spaces/pores of the porous materials. Naturally, the guest gas and water saturations, interfacial tensions and contact angle between the hydrate and liquid phase, size and shape of the irregular internal pores, etc. stand as the authoritative elements. Besides, the salt ions, a key constituent of the saline water, interfere in the formulation of the strong hydrogen-bonded water molecules network that leads to form the hydrate cavities. Albeit they do not participate in the phase transformation, their presence usually reduces the rate of gas hydrate formation and Results formation and decomposition. The  in which, n gg,H is the molar composition of the guest gas (suffix 'gg') in hydrate phase.
Here, the intrinsic rate constant, k is adopted 22 as a function of temperature (T) and activation energy (ΔE a ): k = k 0 exp(−ΔE a /RT).
An important issue is the renewal of reaction surface during the gas hydrate growth and decay. It is taken into account through the time-dependent effective surface area (A e ) of an unconsolidated porous material as, e 0 in which, C is a tuning parameter (constant) that is negative for hydrate growth and positive for decay, and β 0 the surface area adjustment factor. This surface area of the unconsolidated irregular porous sediment (A) is evaluated as: = A e V Tp (2/3) , in which, the scaling factor (e) is 6.2918 27 and the total volume of the porous material (V Tp ) is a function of the porosity (φ) and void volume of the bed (V void ) as: Next issue is the driving force (Δμ) of the gas hydrate formation, growth and decomposition. It is considered as, where, μ w H and μ w L are the chemical potentials of water in the filled hydrate phase and the liquid phase, respectively, which are evaluated from, The fractional occupancy (θ) of each guest gas forming s-I type hydrate structure is a strong function of Langmuir constant (LC) and the fugacity of the guest gas (f j ). Among the terms in Eq. (5), the reference temperature (T 0 ) and the standard chemical potential of water (μ w 0 ) are fixed 28,29 , and the enthalpy (Δh w L ) and molar volume (ΔV w L ) difference between the empty hydrate lattice and liquid water, and the activity of water (a w ) (last component of Eq. (5)) are the function of operating temperature and well pressure 25 . Further, γ w and x w indicate the activity coefficient and concentration of water, respectively.
As far as the second term of last component in Eq. (5) is concerned (porous media contribution), it accounts for the gas hydrate formation in the nanometer sized internal and interstitial pores present in natural porous material (i.e., silica sand), and the associated influence of the surface tension. This term comprises of the radius of a hydrate core present in the internal pores (r core ), the interfacial energy between the planar interfaces (σ ∞ ), the thickness of an interfacial region between the solid and liquid phase (δ) and a linear function of pore radius (k) as: k = a+br pore . Further, the uniform irregularities in the shape of the internal pores are considered through the fractal theory 28 with D f as a fractal dimension 28 . www.nature.com/scientificreports www.nature.com/scientificreports/ The third term of last component in Eq. (5) accounts for the crucial element of the marine hydrate reservoir (i.e., seawater). In which, MW denotes the molecular weight of water, m the concentration of solute species l (i.e., cation, anion or neutral) and φ the osmotic coefficient, which is evaluated from the Pitzer model 30  ). Therefore, these are interrelated quantities, and by knowing any two of them, one can find the third one. For this, we use the following expression as, Solving the differential Eq. (1), we get the respective expression for gas hydrate formation and dissociation as, gg,H In which, α is an adjustable parameter and it is defined as the ratio of the highest value of the net amount of guest gas exhausted during the hydrate formation and the total amount of that gas ideally locked in all cages.
Gas exchange. fundamental mechanism. The difference, in terms of: (i) heat of formation of the CO 2 hydrate and heat of dissociation of the CH 4 hydrate, and (ii) operating pressure and temperature of the formation of CO 2 and CH 4 hydrates, drives the CH 4 -CO 2 (pure/mixed) gas exchange process for the natural gas recovery. This leads to perform double duty i.e., extracting the energy in the form of natural gas and sequestrating the major greenhouse gas contributor (i.e., CO 2 ) in the NGH reservoir. Basically, the pure CO 2 or mixture of the CO 2 -N 2 /H 2 is injected in the existing natural gas hydrate bearing sediments to disturb their equilibrium. Once that gets derailed, the CH 4 hydrate decomposes to the free CH 4 gas and water. At that time, the CO 2 present in the gas phase makes use of these free water molecules to form the CO 2 hydrates 6,9,11 . Thanks to the memory effect of the dissociated water, owing to which, the CO 2 experiences relatively small resistance to form hydrate cages 6,9,11,31 . However, monopoly of the CO 2 (pure/mixed) in forming the hydrates continues until the CH 4 in the gas phase reaches the critical methane concentration (CMC). Afterward, the CH 4 again participates in the hydrate formation 6,9,11,32 . Modeling. To model this complex mechanism, as portrayed in Fig. 2, one needs to consider a twofold role of CH 4 involved in the dissociation and reformation of the gas hydrates. Herein, as a preliminary approach, we are formulating the CH 4 displacement rate during the gas exchange process as, www.nature.com/scientificreports www.nature.com/scientificreports/ Accordingly, the mathematical formulation is proposed as, where, suffixes f and d are used for the formation and decomposition of hydrate, respectively. RG is the replacing gas agent, such as pure CO 2 , mixture of CO 2 -N 2 , CO 2 -H 2 etc. Now, γ, the ratio of the sum of fractional occupancy of the CH 4 in small and large hydrate cages to that of the replacing gas is estimated as, in which, these fractional occupancies are evaluated from the Langmuir type expression, which is detailed elsewhere 28 .
With this concept, the replacement model is obtained as follows, All the specific terms involved in the above equation are elaborated later (see 'Methods' section). This is the final representation of the proposed model to predict the transient behavior of CH 4 -CO 2 (pure/mixed) replacement with porous media in the presence of salt water.
Determining optimal model parameters. The proposed formulation includes total four parameters, namely α, β 0 , C and k 0 . They are inherently linked with the real-time occupancy of guest gases, involved surfaces of porous sediments, renewal in reaction area and rate constant, respectively. These parameters are optimized with the following methodology. optimization problem.
U i n dp 1 gg gg gg Subject to constraints: As shown, the objective function is concerned with the absolute average relative deviation (AARD) of the proposed model predictions (pred) from the actual/experimental (exp) data. U combines a set of optimization parameters as, The objective function is minimized with respect to the proposed model parameters α, β 0 , k 0 and C. For this, we use a nonlinear programming method, namely generalized reduced gradient (GRG). In this method, the nonlinear objective function is linearized at a local solution through the Taylor series expansion. The optimal solution reaches when the partial derivative of an objective function with respect to decision variables is equal to zero 33 . The GRG method is only one of a class of algorithms based on implicit variable elimination. This process begins with the evaluation of reduced gradient vector ∇  f as, www.nature.com/scientificreports www.nature.com/scientificreports/ where, t is the iteration, J and C the constraint gradient sub-matrix of basic (x) and non-basic variables (x ).
, then process terminates, otherwise one has to set descent direction (dd) as: with respect to scalar parameter ω. This nonlinear optimization is proposed to perform effortlessly in the solver tool of Microsoft Excel 2010.
Case study: 2011-2012 Iġnik Sikumi field test. The prediction potential of the proposed formulation will be analysed by evaluating its performance to predict the actual field scale data. For this, we consider the Iġnik Sikumi field test conducted in Alaska North Slope during 2011 and 2012 34 . The primary objective of this project was to investigate the feasibility and real-time performance of the CH 4 -CO 2 /N 2 gas exchange dynamics into the actual hydrate bearing reservoirs. In this test, the injection and production phases are performed through a single vertical well from the temporary ice-pad situated in Prudhoe Bay Unit L-pad. This scheme is usually referred to "huff and puff " method. The process comprises of two principal operational phases: (i) injecting/introducing highly recognized and effective replacing gas mixture of CO 2 (23 mol%) and N 2 (77 mol%) into the reservoir, and (ii) producing CH 4 and recovering CO 2 and N 2 with and without any assistance from the same reservoir. This process was started with the injection of the aforesaid replacing agent at a pressure and temperature of 9.8 MPa and 278.15 K, respectively into the reservoir for 14 days, which is followed by 2.5 days of shut-in soak period. During these 2.5 days, the bottom-hole pressure is reduced from 9.8 to 8.27 MPa. Then, the production/flowback phase was conducted in four stages as, • Stage 1 -Unassisted flowback for 1.5 days • Stage 2 -Jet-pump assisted flowback with P production > P MH stability for 9 days • Stage 3 -Jet-pump assisted flowback with ≅ P P production MH stability for 2.5 days • Stage 4 -Jet-pump assisted flowback with P production < P MH stability for 18.5 days Basically, this production operation was carried out into two broad parts i.e., unassisted (1.5 days) and jet-pump assisted (30 days) flowback with the proper shut-in period for repairing the associated equipment. Further, the jet-pump assisted phase is divided into three stages by manipulating the production/well pressure (P production ) with reference to the methane hydrate stability pressure (P MH stability ). The detailed experimental observations associated with this field test are open to public 35 . performing simulation experiment. We aim to simulate the proposed formulation to interpret the production behavior of the guest gases (CH 4 , CO 2 and N 2 ) involved in the operation of all the stages. For this, we have developed a sequential computer-assisted algorithm as outlined in Fig. 3.
At first, all the components of the proposed formulation need to be analysed before their combined evaluation. They mainly include the total reaction surface area of the involved porous sediment and the driving force of the operation.
evaluating reaction surface area. The well-log analysis made based on the well log and mud log data confirmed 34 that the gas hydrates are present in close proximity to the Prudhoe Bay unit L-pad. Further, the geophysical reservoir characterization of the L-106 well confirmed the presence of four different hydrate bearing sandstone sediments namely, C-2, C-1, D and E. Among these, the homogeneous, thick and highly CH 4 hydrate saturated sediment (i.e., "C-1 sand") has been selected as a primary testing zone. This layer is located in the range of 683.67-692.81 m (2,243-2,273 ft) below the sea level and has 40% porosity, 72% average gas hydrate saturation (S H ) and 28% water saturation (S w ). This region typically has the pressure and temperature of 6.9 MPa and 278.15 K, respectively with the formation breakdown pressure of 10 MPa. For this hydrate bearing sediment, the surface area of the unconsolidated irregular porous sediment is estimated from Eq. (2). In that equation, the main challenge is to estimate the total volume of the porous material (V Tp ) available in this testing zone. For this, we use the value of the void volume and porosity of the bed. Among these, the porosity of the concerned sediment is available 34 . Thus, the accuracy of the evaluated reaction surface is exclusively dependent on the porosity of the bed (here, 40%). Furthermore, the void volume of the bed is projected from the average gas hydrate saturation (here, 72%) and the total volume of methane hydrates (V T,MH ) across the bed as: V void =V T,MH /S H , in which, V T,MH is calculated from the total volume of the recovered methane, assuming that almost all the methane hydrates present in the testing zone is recoverable.
Assessing driving force. The next important component of the proposed formulation is the driving force.
Here, we use the difference in chemical potentials of water in hydrate and liquid phase as a driving force, which is defined in Eq. (3). There are two major components involved, namely μ w H and μ w L [Eqs. (4) and (5)]. Among these, μ w H is a function of the fractional occupancy of the involved guest gas. This fractional occupancy further depends on the fugacity of the guest gas and the Langmuir constant. Here, the temperature dependent Langmuir constant is calculated from the Kihara type potential parameters and the detailed procedure is available elsewhere 22 . Here, the Kihara potential parameters for CH 4 , CO 2 and N 2 include: (i) a (Å) equals 0.3834, 0.6805 and 0.6805, respectively, (ii) σ (Å) equals 3.14393, 2.97638 and 3.13512, respectively, and (iii) ε/K (K) equals 155.593, 175.405 and 127.426, respectively 36 . Now, in the course of four production phases, the well pressure fluctuates, and consequently, the fugacity of each component in the produced gas mixture gets altered, which then affects the hydrate cage occupancy of the involved guest gases -CH 4 , CO 2 and N 2 . Thus, the hydrate phase chemical potential of water changes during the flowback phases. For this, we have evaluated the fugacities of all those three guest gases from the Soave-Redlich-Kwong (SRK) equation of state 37 at each time step using the indigenously developed MATLAB code.
The chemical potential of water in the liquid phase comprises of Δh w L , ΔV w L and a w , which are the function of operating temperature and well pressure. Like the fractional occupancy, these components are evaluated at each time step for the fluctuating well pressure and at the operating temperature. Importantly, a w is computed by accounting the practical issues associated with the two principal constituents (porous sand sediment and seawater) of the marine hydrate reservoirs. As far as the current case study is concerned, the hydrate bearing unit has 72% CH 4 hydrates and 28% water (9% free and 19% bound), leaving no room for the free gases. Naturally, a very less amount of water is available to associate with the injected guest gases to form the new hydrates, and there will be no alternative for the injected CO 2 -N 2 (23-77 mol%) other than following the direct gas exchange mechanism. Besides, from the well log data, the quartz is perceived as a dominant reservoir mineral, and there was no significant mineral dissolution or precipitation of the reservoir minerals into the liquid phase (water) during the operation 16 . Thus, for this case, it is not reasonable to consider their contribution during the simulation run. For simulating the present case reservoir data, the influence of associated water-guest gases and porous material is taken into account. estimating model parameters. Once the calculation of A and Δμ is over, we proceed to estimate the model parameters i.e., α, β 0 , k 0 and C, as mentioned in Fig. 3. We start with the guess value of the tuning parameters and simulate the proposed model to predict the moles of CH 4 , CO 2 and N 2 gases recovery during the production stages. Afterward, the deviation of the predicted response from the field data is used to find the absolute average relative deviation (AARD). We set this AARD as an objective function for the nonlinear optimization problem subject to the constraints stated before. As mentioned, for this, we have used the generalized reduced gradient nonlinear optimization method in the Solver tool inbuilt in Microsoft Excel 2010 setup. This robust optimization technique takes a few seconds to optimize the extensive production data of 31.5 days (unassisted and jet-pump assisted flowbacks). Values of the optimal model parameters are placed in Table 1. The proposed formulation can easily be used for online tuning at diverge geographical conditions. Predicting field data. In this section, the proposed formulation is employed to predict the Iġnik Sikumi field data for natural gas recovery with carbon dioxide sequestration. For this, Figs. 4-6 are produced to investigate the model performance in terms of the dynamic production profile of CH 4 , CO 2 and N 2 , respectively during the four production phases. In this operation 34 , total 215.9 Mscf (6,113.6 m 3 ) of the mixture of replacing agents (CO 2 -N 2 ) was introduced into 9.14 m (30 ft) thick C-1 sand of L-106 bore well at the injection pressure and temperature of 9.8 MPa and close to 278. 15    www.nature.com/scientificreports www.nature.com/scientificreports/ At the start of the production phase, for the unassisted flowback (stage 1-1.5 days), the CH 4 gas recovery rate is slow, and it has dominated the composition in the produced gas mixture. In the second stage (~9 days), the composition of CH 4 in the gas mixture further increases and it reaches around 80 mol% of the total gas produced during this period. Then, for the third (~2.5 days) and fourth stages (~18.5 days), the composition of CH 4 goes on swelling (up to 95 mol%) in the produced gas mixture, while CO 2 and N 2 share a very little contribution. In fact, CO 2 has never exceeded 2 mol% of the total produced stream of the gas mixture 34 . For predicting these outcomes, we use the model having a set of optimized parameters provided in Table 1. With this, the proposed formulation predicts the CH 4 , CO 2 and N 2 gases recovery with only 2.83, 0.84 and 1.67% AARD, respectively. From the prediction accuracy, it becomes obvious that the proposed model is rigorous and versatile enough.
predicting laboratory scale data. Hydrate formation. To investigate the versatility of the proposed formulation, we further compare the performances of the proposed and existing 22 model for the laboratory scale experimental data of CH 4 hydrate kinetics. Figure 7 depicts the comparative performance with reference to the experimental data 38 of water conversion to CH 4 hydrates in the bed of silica sand with 100% water saturation at operating pressure and temperature of 5.75 MPa and 274.5 K, respectively. This setup (Fig. S1) is presented in the supplementary information file and optimal model parameters are listed in Table S1. It is evident that the proposed model performs better than the existing model 22 , which is also obvious from the respective AARD (%) values reported in Fig. 7. This is because the proposed model takes into account many practical aspects involved in hydrate dynamics.  www.nature.com/scientificreports www.nature.com/scientificreports/ Further, we evaluate the proposed model for the CH 4 hydrate growth in the mixture of silica sand and clay. Figure 8(a,b) compare the hydrate growth in the bed of 50% silica sand and 50% clay with 75% and 50% water saturation, respectively. In both the cases, the operating pressure and temperature are 6.10 MPa and 274.5 K, respectively. It is seen in Fig. 8 that the proposed model predicts the CH 4 hydrate growth better than the existing  Hydrate dissociation. Figure 9(a-c) evaluate the performance of the proposed framework comparing with that of the existing model 39 with reference to the experimental data 40 of CH 4 hydrate dissociation in the presence of silica sand, along with pure water, and 1.5 wt% and 3.0 wt% salt solution, respectively. In all these three cases, the operating pressure and temperature are 4.8 MPa and 297.2 K, respectively. The experimental setup is briefly presented in the supplementary file (Fig. S2). The optimal model parameters are documented in Table S2. For all the cases, the proposed model outperforms the existing one 39 , which is quantified through the AARD values reported in Fig. 9. This is because the proposed model has addressed several practical issues related to the porous media and saline water.

conclusion
This work proposes a theoretical formulation to portray the gas hydrate formation, decomposition and replacement dynamics in the actual hydrate reservoirs. This rigorous model formulates the two-fold CH 4 -CO 2 (pure/ mixed) exchange mechanism, addressing several practical aspects specific to the basic elements (porous media and saline water) of the hydrate reservoirs. Potential of the proposed framework is analysed by predicting the amount of CH 4 produced, and CO 2 and N 2 recovered in the Iġnik Sikumi field test conducted on the Alaska North Slope in 2011 and 2012. This field test starts with the injection of CO 2 -N 2 (23-77 mol%) gas mixture into the well for 14 days, which is then followed by 31.5 days of CH 4 production and CO 2 -N 2 recovery from the same bore well. For this, we perform the simulation experiments with the conditions at which the field scale experiment was exactly conducted. Thanks to the nonlinear optimization technique, using which, we have identified the optimized parameter set online for all the production phases. With this, the model has shown a promising performance in predicting the field scale data. It is also quantified based on the percent AARD obtained for CH 4 , CO 2 and N 2 as 2.83, 0.84 and 1.67, respectively with reference to the Iġnik Sikumi field data. To show the versatility of the proposed framework, further it is successfully tested with the laboratory scale data for CH 4 hydrate formation and decomposition in the presence of pure and salt water, and porous media. This approves the online use of the proposed theory to interpret the actual field scale exploration of the naturally occurring hydrates through the replacement technique.

Methods
Modeling. The amount of water associated with the hydrate and liquid phase is expressed in terms of n CH4,H , η H and n RG,H . With this, Eq. (10) is rewritten as, The above equation is simplified to, CH ,H CH , H 4 4 in which, the time dependent parameters P and Q have the following representation, Since at t = 0, n gg,H = n 0 , the above equation yields,