Development and Application of a Simulator for Simulating the Behaviors of a Geological System When Replacing CH₄ from Hydrate-Bearing Reservoirs by CO₂

Abstract

Conventional techniques for hydrate production may cause the deconstruction of hydrate, changing the geomechanical stresses of the reservoir, which could trigger the subsidence of the seafloor. A new method for replacing CH₄ from the hydrate lattice by CO₂, without damaging the mechanical structure of sediment, has been proposed. This approach can achieve both the objectives of long-term CO₂ sequestration and the safe production of CH₄ from hydrates. By coupling the Chen-Guo model into Tough+Hydrate V1.5, an updated simulator CO₂-EGHRSim V.10 (CO₂ Enhanced Gas Hydrate Recovery simulator) was developed in this work to describe the replacing processes of CH₄ from the hydrate lattice by CO₂ and to evaluate the storage potential of CO₂ and the recovery efficiency of CH₄ from the hydrate-bearing reservoirs. The developed simulator was verified using measured data obtained from laboratory experiments. The verification suggested that CO₂-EGHRSim performed well in predicting the replacing processes of CH₄ with CO₂. The simulator was applied to calculate the CO₂ storage potential combined with the CH₄ recovery from hydrates at the site of Iġnik Sikumi on the North Slope of Alaska. The simulated results indicated that the CO₂–CH₄ exchange mostly occurred inside the gas plume, and the CO₂ hydrate was only present around the production well. The simulated CO₂ storage ratio was 0.58, and the CH₄ recovery efficiency was 25.95%.

1. Introduction

Carbon geological utilization and storage (CGUS) is a promising technological option for CO₂ reduction. CO₂-enhanced gas hydrate recovery (CO₂-EGHR) is a technique of CGUS for CO₂ geological utilization. The technique allows CH₄ to be safely generated from the hydrate-bearing reservoir and simultaneously allows CO₂ to be trapped by the geological formation in the form of hydrate [1].

CO₂-EGHR consists in injecting the pure CO₂ or a mixture composed of CO₂ and other gases into natural gas hydrate reservoirs to extract CH₄ from the hydrate lattice [2,3,4]. The feasibility of the replacing process is owed to the fact that at temperatures and pressures typical for geological CH₄ hydrate formation, a CO₂ hydrate exhibits more thermodynamic stability than a CH₄ hydrate [5]. Laboratory experiments have validated the thermodynamic feasibility of the CO₂–CH₄ exchanging method [6]. Due to the lack of widespread approval, this method has not been implemented on a field scale. In addition, carrying out the CO₂-EGHR at the field scale is too expensive to promote its wide application. Reservoir simulations are recommended to estimate the potentials both of CO₂ sequestrated in the solid hydrate phase and of CH₄ recovered from the natural gas hydrates [7].

The simulators used to model the complicated behavior of CO₂ replacing CH₄ from the hydrate deposits have rarely been reported in previous works. A previous work [8] proposed a CO₂ gas injection numerical simulator based on the existing numerical simulator STOMP-HYD [9] and was applied to the Iġnik Sikumi field trial [10]. A simulator called STARS was used for CH₄/CO₂/N₂ guest molecule exchange replacement in the presence of pure water [11]. The previous simulated results [7] implied that the CH₄ − (N₂ + CO₂) exchanging process was complicated and was strongly sensitive to the environmental conditions (e.g., pressure, temperature, salinity, etc.).

Mostly, for a natural gas hydrate, the caged gas in the hydrate lattice is methane; however, a certain percent of other gases could be captured in the cage formed by water molecules [12]. The guest gas cage in the hydrate lattice decides the hydrate structure types, e.g., CH₄ usually forms a type I hydrate, but C₃H₈ easily forms a type II hydrate. CO₂ sometimes forms a type I hydrate and sometimes forms a type II hydrate [13]. Currently, rare simulators are capable of distinguishing the hydrate structure types of hydrates. The complex and susceptible process promotes the demands for a sophisticated numerical simulator to further investigate the replacing processes of CH₄ from a hydrate-bearing reservoir.

To this end, a comprehensive simulator was developed in this work and was used to simulate the replacing behavior of CH₄ from the geological hydrate-bearing reservoir by CO₂. The Chen-Guo thermodynamic model was coupled into the existing Tough+Hydrate V1.5 to estimate the conditions where the mixed-gas hydrates can exist stably. The thermodynamic properties (density, specific enthalpy, specific internal energy, and viscosity) of the liquid phase, gas phase and hydrate phase were estimated using empirical models. The simulator was validated using the data of experiments carried out in the lab and was then applied to a field trial test to estimate the potentials of CH₄ recovery and of CO₂ storage. By discussing the simulated results, main conclusions were drawn.

2. Methodologies and Simulator Development

To perform simulations of extracting CH₄ from the hydrate-bearing reservoir by injecting pure CO₂ or a mixture composed of CO₂ and other gases, a robust simulator is necessary. Tough+Hydrate V1.5 is an open-source simulation tool capable of modeling the complicated behaviors of a hydrate-bearing geosystem [14]. However, only pure CH₄ hydrate can be considered in the released version of Tough+Hydrate 1.5. In this study, under the framework of Tough+Hydrate V1.5, a simulator, CO₂-EGHRSim V1.0, was developed to characterize the behaviors of geological systems bearing hydrates formed by pure CH₄, CO₂, N₂, C₂H₆ and C₃H₈, or their mixtures.

2.1. Mixed-Gas Hydrate Model

Two types of models for predicting the formation/dissociation of a gas hydrate are commonly used, i.e., the kinetic model and the thermodynamic equilibrium model. A previous work [15] concluded that the kinetic model played an important role in predicting the gas production rate at a reactor core; however, the latter one is widely employed to forecast the production of gas on a reservoir scale that is controlled by mass and heat transfer. Chen-Guo’s model includes these two types to calculate the phase boundary and the hydrate formation rate [12]. In this study, the reservoir scale was the focus; therefore, the thermodynamic type of Chen-Guo’s model was used to forecast the conditions (temperature and pressure) in which a mixed-gas hydrate forms. The Chen-Guo model in general is derived from the prediction theory of hydrate stability, which is based on equating the fugacity of gas species i in the gas phase and in the hydrate phase:

$$f_i=x_i^{*}{f}_i^0{\left(1-\sum _j{\theta }_j\right)}^{\alpha }$$

(1)

$$\sum _i{x}_i^{*}=1.0$$

(2)

where $f_i$ is the fugacity of gas species i in the gas phase which forms the basic hydrate i, and in this study, it was computed using SRK-EOS [16]; $x_i^{*}$ denotes the mole fraction of the i basic hydrate, which is formed by the gas i; $\alpha$ is a hydrate-structure-dependent constant, where for structure I hydrate, $\mathsf{\alpha }=1/3$ and for structure II, $\mathsf{\alpha }=2$; $f_i^0$ labels the fugacity of the formed hydrate i equilibrating with a pure basic hydrate and is a multifactor sensitive parameter:

$$f_i^0=f_i^0\left(T\right)·f_i^0\left(P,T\right)·f_i^0\left({\alpha }_W\right)$$

(3)

where ${\alpha }_W$ is the activity of water; $f_i^0\left(T\right)$ is a function of T, and can be expressed as:

$$f_i^0\left(T\right)=\exp \left[\frac{-\sum _j{A}_{ij}{\theta }_j}{T}\right]·\left[A_i^{\prime }\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{B_i^{\prime }}{T-C_i^{\prime }}\right)\right]$$

(4)

where $A_{ij}$ are the binary interaction parameters indicating the relations of the guest gases encaged by the small and large cavities. The values of $A_{ij}$ are listed in Table 3 of [12]. A’, B’ and C’ are Antoine’s constants for different guest gases and structure types and are summarized in Table 3-18 and Table 3-19 of the work [12].

In Equation (3), $f_i^0\left(P,T\right)$ can be calculated using:

$$f_i^0\left(P,T\right)=\exp \left[\frac{\beta P}{T}\right]$$

(5)

where $\beta$ is a hydrate-structure-dependent constant: $\beta =4.242 K {MPa}^{-1}$ for a structure I hydrate and $\beta =10.224 K {MPa}^{-1}$ for a structure II hydrate.

In Equation (3), $f_i^0\left({\alpha }_W\right)$ is given by:

$$f_i^0\left({\alpha }_W\right)={\alpha }_W^{-1/{\lambda }_2}$$

(6)

where λ₂ is a hydrate-structure-dependent parameter, 3/23 for structure I, and 1/17 for structure II, respectively. In this study, the activity of water was assumed as unity for all the system temperatures.

In Equations (1) and (4), θ_j is the fraction of the small cavities that the gas species j occupy and is obtained through the Langmuir adsorption theory:

$${\theta }_j=\frac{f_j{C}_j}{1+\sum _k{f}_k{C}_k}$$

(7)

where the Langmuir constant C for the guest gas can be obtained in [12].

2.2. Estimation of Thermodynamic Properties for Phases

In a hydrate-bearing system, four phases are defined, i.e., the gas phase, aqueous phase, hydrate phase and ice phase. The main thermodynamic parameters (such as density, enthalpy, viscosity, etc.) of a mixed-gas hydrate-bearing system are different with those of a pure CH₄ hydrate system (as can be found in [14]) and are discussed in the following sections.

2.2.1. Thermodynamic Properties of the Gas Phase

Density of the Gas Phase

To estimate the gas phase density, the following expression was employed:

$${\rho }_{mix}=\frac{P·{MW}_{mix}}{Z_V·R·T}$$

(8)

where ${\rho }_{mix}$ is the density of the gas phase, in kg∙m⁻³; the molecular weight of the gas phase is ${MW}_{mix}=\sum _i{y}_i·{MW}_i$, where MW_i is the molecular weight of individual gas species i; and Z_V is the compressibility factor of the vapor phase calculated using SRK EOS.

Enthalpy of the Gas Phase

The enthalpy of the gas phase is determined using the equation of state by combing the departure function for the mixture with the enthalpy of an ideal gas mixture. The enthalpy of a gas mixture is evaluated according to the volume fractions of gases in the mixture. A model (Equation (9) from [17]) was used to calculate the specific enthalpy of the gas mixture.

$$h_{mix}^{real}=h_{mix}^{id}+h_{mix}^{Dep}+\sum y_i{h}_i^{ref}$$

(9)

where $h_{mix}^{id}$ is the ideal gas enthalpy of the gas mixture in J kg⁻¹, which can be expressed by:

$$h_{mix}^{id}=\sum y_i{h}_i^{id}$$

(10)

In Equation (10), based on the specified reference temperature T_ref, the specific enthalpy of gas i at the desired temperature T is calculated by:

$$h_i^{id}=R\sum _{j=1}^5\frac{A_{i,j}\left(T^j-T_{ref}^j\right)}{j}$$

(11)

where R = 8314.56 J K⁻¹ kmol⁻¹. A_i,j are regressed constants from experimental data (

Table 1. Regressed constants of Equation (11).

No.	Gas	Ai,1	Ai,2 (×10−3)	Ai,3 (×10−5)	Ai,4 (×10−8)	Ai,5 (×10−8)
1	CH4	4.568	−8.975	3.631	−3.407	1.091
2	CO2	3.259	1.356	1.502	−2.374	1.056
3	N2	3.539	−0.261	0.007	0.157	−0.099
4	C2H6	4.178	−4.427	5.66	−6.651	2.487
5	C3H8	3.847	5.131	6.011	−7.893	3.079

).

In Equation (9), $h_{mix}^{Dep}$ is the departure enthalpy, estimated using SRK-EOS [16] in this work.

$$h_{mix}^{Dep}=RT\left(Z-1\right)-\frac{a_{mix}}{b_{mix}}\left(1-\frac{T{a}_{mix}^{\prime }}{a_{mix}}\right)\mathrm{l}\mathrm{n}\frac{V+b_{mix}}{V}$$

(12)

where $a_{mix}^{\prime }=\raisebox{1ex}{$d{a}_{mix}$}\!\left/ \!\raisebox{-1ex}{$dT$}\right.$ can be obtained by:

$$\left(\frac{d{a}_{mix}}{dT}\right)=\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{y}_i{y}_j\left(1-k_{ij}\right)\left(\sqrt{\frac{a_j}{a_i}}\left(\frac{d{a}_i}{dT}\right)+\sqrt{\frac{a_i}{a_j}}\left(\frac{d{a}_j}{dT}\right)\right)$$

(13)

$$\frac{d{a}_i}{dT}=\frac{-m_i{a}_i}{\sqrt{{\alpha }_{i}T{T}_{ci}}}$$

(14)

In Equation (9), y_i is the volume fractions of gas species i in the mixture excepting water vapor. $h_i^{ref}$ is the enthalpy of gas species i in its reference phase from the elemental species at P_ref and T_ref, which are generally taken as 1.0 bar and 273.15 K, respectively. The reference phase used in these calculations is the ideal gas.

Due to the polar property of water, in SRK-EOS, the enthalpy of steam is not included when calculating the specific enthalpy of the mixture. Alternatively, the steam enthalpy is added to the specific enthalpy of the real gas:

$$h_{Gas}=\left(1.0-x_{steam}\right)h_{mix}^{real}+x_{steam}{h}_{steam}$$

(15)

where x_steam is the mass fraction of steam in the gas phase, and h_steam is the specific enthalpy of steam and is estimated using a more accurate model, SUPST, proposed by [18].

Viscosity of the Gas Phase

The viscosity is another important property of the gas phase and is calculated using a function at a given pressure, temperature, compressibility, and composition according to the friction theory developed by [19] based on the method of [20]. The mixture viscosity is given by:

$${\mu }_{mix}={\mu }_{0,mix}+{\mu }_{f,mix}=\exp \left[\sum _{i=1}^N{y}_i\mathrm{l}\mathrm{n}\left({\mu }_{0,i}\right)\right]+{\mu }_{f,mix}$$

(16)

where y_i is the mole fraction of gas i in the mixture, and ${\mu }_{0,i}$ is given by:

$${\mu }_{0,i}=40.785\frac{\sqrt{M_i·T}}{V_{c,i}^{2/3}·{\Omega }^{*}}·F_{c,i}$$

(17)

The parameters involved in Equations (16) and (17) are referred to in [19].

2.2.2. Thermodynamic Properties of the Hydrate Phase

In the new package CO₂-EGHRSim V1.0, five gas species are currently accommodated. The thermodynamic parameters of mixed-gas hydrates are significantly different from those of a pure CH₄ hydrate.

Density of the Hydrate Phase

In this study, the typical gases forming hydrates included CH₄, CO₂, N₂, C₂H₆ and C₃H₈. The density of the hydrate formed by the typical pure gases can be found in a previous work [12] and are listed in

Table 2. The density of hydrates formed by typical pure gases (273.15 K).

Gas	CH4	CO2	N2	C2H6	C3H8
Molecular weight (g/mol)	16.04	44.01	28.04	30.07	44.09
Density of hydrate (kg/m3)	910.0	1117.0	995.0	959.0	866.0

.

For the hydrate formed by a gas mixture, the density can be obtained by summing the volumes of individual basic hydrates:

$$\frac{1}{{\rho }_{MGH}}=\sum _{i=1}^N\frac{x_{i,GBH}}{{\rho }_{i,PGH}}$$

(18)

where ${\rho }_{MGH}$ is the density of the hydrate formed by the gas mixture, in unit of kg m⁻³; N is the number of gases forming hydrates; $x_{i,GBH}$ is the mass fraction of basic hydrate i and ${\rho }_{i,PGH}$ is the corresponding density of basic hydrate i (Table 1), in kg m⁻³.

In Equation (18), the mass fraction of basic hydrate i was estimated using:

$$x_{i,GBH}=\frac{x_i^{*}\left({MW}_{i,G}+N_{i,hn}^{type}{MW}_{H2O}\right)}{\sum _{j=1}^N{x}_j^{*}\left({MW}_{j,G}+N_{i,hn}^{type}{MW}_{H2O}\right)}$$

(19)

where $x_i^{*}$ is the mole fraction of basic hydrate i obtained in Equation (2); MW_i,G is the molecular weight of gas i (Table 2), and MW_H2O is molecular weight of water; and $N_{i,hn}^{type}$ is the hydration number of basic hydrate I formed by gas i. For different structures, $N_{i,hn}^{type}$ are generally estimated using:

$$\begin{gathered} N_{i,hn}^I=\frac{23}{3+{\theta }_{i,1}} \\ {N}_{i,hn}^{II}=\frac{17}{1+2{\theta }_{i,1}} \end{gathered}$$

(20)

$N_{i,hn}^I$ and $N_{i,hn}^{II}$ are the hydration numbers of structure I and II hydrates; θ_i,1 is the fraction of the small cavities occupied by gas species i and is derived from Equation (7).

Enthalpy of the Hydrate Phase

During the behaviors of the hydrate-bearing geosystem, the specific enthalpy is used to calculate heat transfer. Regarding the hydrate phase, the specific enthalpy is estimated through the following expression:

$$\begin{gathered} h_{mix}=\sum _{i=1}^N{x}_i^{*}{h}_i \\ {h}_i=\underset{T_{ref}}{\overset{T}{\int }}{C}_{i}dT \end{gathered}$$

(21)

where $x_i^{*}$ is the mole fraction of basic hydrate i in the mixed gas hydrate; h_i is specific enthalpy of basic hydrate i; and T_ref is defined as the reference temperature, 273.15 K; C_i is the heat capacity of basic hydrate i, in J kg⁻¹ K⁻¹. Regarding the involved pure gas hydrate in this study, the heat capacities of pure CH₄ hydrate (2098.00), C₂H₆ hydrate (2064.00), C₃H₈ hydrate (2075.00), CO₂ hydrate (1698.00) and N₂ hydrate (1543.00) were obtained from the work of [21,22].

Dissociation Heat of the Hydrate

A uniformed model [23] was adopted to calculate the dissociation heat of the mixed-gas hydrate:

$$∆H_D=a+bT$$

(22)

where ∆H_D denotes the heats of dissociation for hydrates, in Cal∙gmol⁻¹; a and b are hydrate-specific coefficients regressed from experimental data, listed in

Table 3. Coefficients in Equation (22) for different gas hydrates [23].

Gas	a × 10−3	b	T Range (K)
CH4	6.534	−11.97	248.15–273.15
	13.521	−4.02	273.15–298.15
CO2	9.29	−12.93	248.15–273.15
	19.199	−14.95	273.15–284.15
N2	4.934	−9.04	248.15–273.15
	6.188	18.37	273.15–298.15
C2H6	8.458	−9.59	248.15–273.15
	13.254	−15.00	273.15–287.15
C3H8	7.06	−4.09	248.15–273.15
	−37.752	250.09	273.15–278.15

; T is the temperature in K.

3. Code Verification

The accuracy and reliability of CO₂-EGHRSim V1.0 should be verified before it can be applied to a field site. Based on the data from laboratory experiments, the developed simulator was verified. The precision of the simulation results was determined using the average absolute relative deviation (AARD, Equation (23)).

$$\mathrm{A}\mathrm{A}\mathrm{R}\mathrm{D}=\frac{1}{N}\sum _{i=1}^N\left(\frac{\left|{Data}_i^{Exp}-{Data}_i^{Sim}\right|}{{Data}_i^{Exp}}\right)$$

(23)

where N is the number of experimental data points, ${Data}_i^{Exp}$ is the measured data, and ${Data}_i^{Sim}$ is the simulated data.

3.1. Verification of the Prediction of Hydrate Formation Conditions

In this section, the predicted formation conditions and structure types of the mixed-gas hydrates using the developed simulator were compared with the values obtained from the laboratory experiments. The experimental measurements of hydrate conditions for pure gases, binary, ternary, and quaternary mixtures were conducted by [24].

Under high-pressure conditions, the experiments were conducted in an autoclave with a capacity of approximately 500 mL. The container was made of corrosion- and acid-resistant stainless steel, which was designed to be able to withstand pressures of up to 25 MPa. Pure CH₄, N₂, CO₂, C₂H₆, C₃H₈ and their arbitrary mixtures were used to form hydrates. Nineteen cases, including four pure cases, five binary cases, five ternary cases and five natural-gas cases, were designed in the laboratory experiments. Only thirteen cases (

Table 4. Cases used in the verifications.

Case	Composition (Mol%)
	CH4	CO2	N2	C2H6	C3H8
P1	100	/	/	/	/
P2	/	/	100	/	/
P3	/	/	/	100	/
P4	/	/	/	/	100
B1	89.26	/	10.74	/	/
B2	90.47	/	/	9.53	/
B3	97.07	/	/	/	2.93
B4	/	/	/	85.15	14.85
T1	94.97	5.00	0.03	/	/
T2	90.93	/	4.18	4.89	/
T3	84.52	/	/	12.55	2.93
Q1	89.40	8.09	0.02	/	2.49
Q2	89.60	5.25	0.02	5.13	/

Note: P—pure gas case; B—binary gas case; T—ternary gas case; Q—quaternary gas case; /—was not contained.

) were employed to verify the simulator for predicting the formation conditions of hydrates.

The simulated and measured formation pressures of hydrates at given temperatures were compared and are plotted in Figure 1.

Figure 1 shows that for both pure-gas cases and for mixed-gas cases, the predicted formation pressures of hydrates at given temperatures using the developed simulator (solid lines with different colors) were in good agreement with the measured data (indicated by varying symbols). The AARD of all case simulations were calculated using Equation (23) and are listed in

Table 5. AARD calculated for all case simulations.

Case	Number of Data Points	AARD (%)
P1	30	1.17
P2	5	1.82
P3	15	1.38
P4	10	4.22
B1	6	2.11
B2	8	4.90
B3	4	4.52
B4	6	0.65
T1	6	1.28
T2	7	4.33
T3	5	4.73
Q1	5	4.94
Q2	6	4.46

. Analyzing the calculated AARD showed that all AARD values obtained in the thirteen case simulations were below 5.0%. The case simulation of Q1 yielded the maximum AARD of 4.94%, and the case simulation of B4 gave the minimum AARD of 0.65%. The AARD derived from the pure, binary and ternary gas cases were lower than those from the quaternary cases.

Both the comparison of simulated values with measured data and the calculated AARD suggested an acceptable accuracy of the developed simulator CO₂-EGHRSim for estimating the hydrate formation conditions.

3.2. Verification on the Replacement Process of CH₄ from Hydrate with CO₂

In this section, to verify the simulator CO₂-EGHRSim, the simulated replacement process of CH₄ from the hydrate phase with CO₂ and N₂ were compared with the data from laboratory experiments [25]. The experiments were carried out in a high-pressure and low-temperature reactor (ϕ 44 × 200 mm) in which silica sand and distilled water were filled, acting as sediments. After the injection of the gas mixture, the closed reactor was put into a thermotank. The first step of the experiments was to form a pure CH₄ hydrate by changing P and T, and then CO₂ and N₂ or their mixture were injected to replace CH₄ from the CH₄ hydrate. The volume fractions of CH₄, CO₂ and N₂ were measured using a gas chromatograph. Two group of experiments were selected in this study to verify the simulator. The first group was carried out at a fixed pressure of 14.0 MPa with a constant injectant composition of 53% CO₂ + 47% N₂, and the temperature of this group changed from 1.0 to 10.0 °C; the second group was performed at a constant temperature of 1.0 °C, the injectant compositions were set to 100% CO₂ + 0.0% N₂, 53% CO₂ + 47% N₂ and 19% CO₂ + 81% N₂ to keep a fixed partial pressure of CO₂ of 3.4 MPa in the system; thus, the second-group experiments were conducted at pressures of 3.4 MPa, 6.5 MPa and 15.8 MPa, respectively. The experiment settings are listed in

Table 6. Cases used in the verification.

Case No.	Upper Block (Gas Phase)				Lower Block (Aqueous + Hydrate, and Saq = 0.41)
	P(MPa)	T (°C)	CO2 (VFr)	N2 (VFr)	P (MPa)	T (°C)	BHMFr
1	14.00	1.00	0.53	0.47	14.00	1.00	CH4: 1.00 CO2: 0.00 N2: 0.00
2		4.00				4.00
3		8.00				8.00
4		10.00				10.0
5	3.40	1.00	1.00	0.00	2.50	1.00
6	6.50		0.53	0.47
7	15.80		0.19	0.81

VFr: volume fraction; MFr: mole fraction; BHMFr: mole fraction of basic hydrate.

. The replacement experiments were carried out for different days. The longest time of the experiments selected in this study for verifying the simulator was 10 days. To observe the behavior of the simulator, we extended the simulation time to 11.5 days. The recovery ratio of CH₄ was employed to estimate the replacement effects and can be calculated using:

$$Re=\frac{n_{CH4,G}}{n_{CH4,Init}}$$

(24)

where Re is the CH₄ recovery ratio, n_CH4,G is the moles of CH₄ in the gas mixture, and n_CH4,Init is the moles of CH₄ initially present in the hydrate phase.

To match the experimental conditions, a two-block conceptual model was built (Figure 2) to characterize the reactor used in the experiments. The diameter of the two blocks was 44.0 mm. The upper block of the conceptual model was a 160.0 mm high block saturated by the gas mixture. The lower block was a 40.0 mm high block with distilled water and silica sand.

In this study, the CH₄ recovery ratios of the two groups of experiments (53% CO₂ + 47% N₂ gas mixture at temperatures of 1.0 °C, 4.0 °C, 8.0 °C and 10.0 °C; CO₂, 19%CO₂ + 81%N₂, and 53% CO₂ + 47% N₂ gas mixtures at a temperature of 1.0 °C and pressures of 3.4 MPa, 6.5 MPa and 15.8 MPa) were simulated using CO₂-EGHRSim V1.0 and compared with the experimental data.

The upper block was designed to represent the headspace of the reactor. The upper block was set to hold a higher porosity of 0.98 and a higher permeability of 1.0 × 10⁻¹⁰ m². In contrast, the lower block acted as the sediment where the formation of a hydrate was expected to occur and possessed a lower porosity of 0.39 and a lower permeability of 1.0 × 10⁻¹³ m². All other geophysical parameters were set according to the experiments conducted by [25] and are listed in Table 6.

After 11.5-day simulations, both the simulated and the measured recovery ratios of CH₄ by CO₂ were obtained and are plotted in Figure 3. It can be found that all the simulated CH₄ recovery ratios (solid lines in different colors) were very close to the measured values (symbols with different colors) in the experiments. At the beginning of the replacements, the measured recovery ratios of CH₄ were slightly lower than the simulated values, which could be explained by the utilization of the thermodynamic model for predicting the formation of hydrates.

As indicated by the experimental and simulation data (in Figure 3), the recovery ratio increased with the temperature (Figure 3a) and CO₂ molar fraction (Figure 3b). A higher temperature increased the dissociation of hydrate releasing more CH₄ than a lower temperature, which raised the recovery ratio of CH₄ from the hydrate phase by CO₂. As for the molar fraction of CO₂, please note that the pressures for cases with different CO₂ molar fractions were different. The pressure for the case with pure CO₂ was 3.4 MPa, the pressure for the case with 53% CO₂ + 47% N₂ was 6.50 MPa, and the pressure for the case with 19% CO₂ + 81% N₂ was 15.80 MPa. Those cases were designed to keep the partial pressure of CO₂ at 3.4 MPa. Thus, the case with a higher molar fraction of CO₂ held a lower system pressure. At the beginning of the replacement, in the case with a higher pressure, the amount of formed CO₂ hydrate was larger, which hindered the dissociation of CH₄ hydrate and lowered the recovery ratio [25].

The AARDs were calculated for all cases based on the measured data and are listed in

Table 7. AARD calculated in all case simulations.

Case No.	Number of Data Points	AARD (%)
1	6	2.61
2	7	3.22
3	8	4.71
4	8	4.84
5	7	3.73
6	6	2.93
7	6	2.61

. It can be concluded that the simulated recovery ratios of cases with different temperatures and gas compositions were in good agreement with the measured data with a maximum AARD of 4.84% lower than 5.00%.

Based on the laboratory experiments, the accuracy of the developed simulator CO₂-EGHRSim V1.0 was validated through verifying the prediction of hydrate formation conditions and the description of the replacement of CH₄ from the hydrate phase by CO₂. The verification results suggested a high accuracy and an acceptable reliability.

4. Application of the Simulator to a Field Trial at the Iġnik Sikumi Site

After verification, the developed CO₂-EGHRSim V1.0 simulator was applied to a field trial at the Iġnik Sikumi site on the North Slope of Alaska. The field test was carried out to investigate the CO₂–CH₄ exchange on a large field scale, and to estimate the storage potential of CO₂ and recovery efficiency of CH₄ from the hydrate-bearing reservoir.

4.1. Brief Description of the Iġnik Sikumi Site

In 2011 and 2012, ConocoPhillips conducted a field trial at the Iġnik Sikumi site in the Prudhoe Bay Unit on the Alaska North Slope (Figure 4) to analyze the feasibility of CH₄ recovery by injecting CO₂ into the hydrate-bearing reservoirs. The test trial site was selected according to the comprehensive analyses of the accessibility, proximity to the North Slope infrastructure and the present possibility of a hydrate-bearing reservoir. The reservoirs mainly made of clastic sandstones of Sagavanirktok formation, which usually possesses a high porosity and permeability, were selected as the target injection layer for the field trial. The injection well perforated over a 9.14 m (30 ft) interval (683.67 to 692.81 m, 2243 to 2273 ft). The field trial proceeded through four stages, including (1) a 14-day injection; (2) 3-day soak operations with pressure recovery; (3) a 1.5-day unassisted flow-back; and (4) a 30-day jet-pump-assisted flow-back.

A controlled mixture of N₂ and CO₂ with mole fractions of 77% and 23% [26] was injected with a constant pressure of 9.79 MPa (1420 psia) at a depth of 695.86 m (2283.0 ft) for 306 h, with vertical pressure gradient of 9817.32 Pa/m (0.434 pisa/ft). The production of the field trial was split into a 7-day low-flow period, a 2.5-day high-flow period, and a 19-day increasing-flow period.

4.2. Conceptual Model and Boundary Conditions Setting

To investigate the exchanging processes of CH₄ with CO₂ from hydrates, in this study, a domain with a vertical length of 15.24 m (from 680.92 m to 696.16 m) was selected. The injection/production well perforated from 683.67 m to 692.81 m. The radial extension of the domain began with the well casing (0.1143 m) and ended at a distance of 609.6 m.

To describe the domain for the CO₂–CH₄ exchange, a 2D radial conceptual model was constructed (Figure 5). The model was 15.24 m thick and 609.6 m long in radial distance. The model was evenly divided into 25 layers in the vertical direction, and logarithmically divided into 50 blocks in the radial direction. Thus, the whole model was discretized into 1250 blocks.

The outer radial boundary, the top boundary and the bottom boundary of the model were considered as the Dirichlet boundary conditions. The temperature and pressure at those boundaries were fixed at the initial values where the CH₄ hydrate was stable. The boundary condition settings are shown in Figure 5.

4.3. Initial Conditions and Geophysical Parameters

A gas mixture with a fixed composition of 23% CO₂ and 77% N₂ was injected from the center of the computational domain. At the perforation, the injection pressure was kept at about 9.86 MPa, and the corresponding injection temperature was kept at about 5.39 °C. A vertical pressure gradient of 9.82 kPa∙m⁻¹ [10] was used to obtain the pressure distribution along the injection well perforation with a constant pressure of 9.86 MPa at a depth of 695.00 m.

In the computational domain, the hydrate and free water initially coexisted. According to the primary variables’ setting, P, T and the saturation of the hydrate were required. These three primary variables varied vertically for each 0.6096 m thick layers. The initial values of the primary variables (P, T and the initial saturation of the hydrate) were taken from the logging data reported in the Table 4 of the Final Technical Report by [10].

The geophysical properties (porosity, permeability, etc.) were taken from the settings in Table 1 of the work by [10]. For the injection well, grid blocks with a porosity of 1.0 and a permeability of 1.0 darcy were used to represent the wellbore [27]. Actually, the porosity and permeability were not suitable for the wellbore. However, in the simulator CO₂-EGHRSim updated from Tough+Hydrate, the multiphase fluid flow was assumed to follow Darcy’s Law. In order to mimic the well, some grid blocks with a higher porosity and permeability were set. Note that those grid blocks just provided injectants (mixtures composed of CO₂, N₂ and CH₄) with a constant injection pressure. At the beginning of the simulations, only the CH₄ hydrate was assumed to be present in the hydrate-bearing reservoir. The initial distribution of the hydrate saturation is plotted in Figure 6 based on the initial hydrate distribution data reported in [10].

4.4. Results and Discussion

The Iġnik Sikumi field trial mainly included an injection stage and a production stage. During the injection stage, the gas mixture of 23% CO₂ and 77% N₂ was injected for about 306 h with a constant injection pressure. The complete production period was about 30 days. The injectivity of CO₂ and the recovery efficiency of CH₄ are evaluated and simultaneously, the CH₄–CO₂ exchange processes in the recovery of CH₄ are investigated in the following sections.

4.4.1. Evaluation of Gas Mixture Injectivity

At the end of the 306 h injection, the amounts of injected CO₂ and N₂ were evaluated and are plotted in Figure 7. The simulated gas mixture injected into the computational domain weighed 0.79 × 10⁴ kg, including 0.25 × 10⁴ kg of CO₂ (green dash line in Figure 7) and 0.54 × 10⁴ kg of N₂ (blue solid line in Figure 7), respectively. In the field trial, the measured data showed that the injected gas mixture weighed 0.80 × 10⁴ kg with CO₂ and N₂ masses of 0.25 × 10⁴ kg (green triangles in Figure 7) and 0.55 × 10⁴ kg (blue rectangles in Figure 7), respectively. The relative errors were 2.07%, 1.04% and 2.53% for the injected gas mixture, the injected CO₂ and the injected N₂, respectively.

4.4.2. Changes of Physical Properties Due to Injection of CO₂ and N₂

The saturation of the gas phase after injection is shown in Figure 8a. After the 306 h migration, the front of the gas phase reached about 8.80 m from the outer area of the injection well. Under suitable conditions, pure- or mixed-gas hydrates formed. The formation of the hydrates was an exothermic reaction, which increased the temperatures of the zones where hydrates formed (Figure 8b). It was easier to form a CO₂ hydrate than a CH₄ hydrate at the same temperature and pressure, and the swapping process of CH₄ with CO₂ was spontaneous [28,29]. Figure 8c gives the distribution of the hydrate saturation within the computational domain at the end of the injection. Compared with the initial state (Figure 6), the distribution of the hydrate saturation obviously changed near the injection well. The saturation of the hydrate phase increased near the injection well due to the formation of CO₂ hydrate and N₂ hydrate. The saturation reached 0.83 near the injection well where the initial saturation was about 0.62. The increment in hydrate saturation raised the temperature of the zone where the formation of hydrate dominated (as shown in Figure 8b).

In the hydrate phase, the mass fractions of basic hydrates formed by individual gases were calculated and are depicted in Figure 8d–f. The mass fraction of CH₄ hydrate declined near the injection well from 1.00 to 0.02 (Figure 8d) due to the replacement by CO₂ hydrate or N₂ hydrate. The mass fraction of CO₂ hydrate increased from 0.00 to 0.25 m near the injection well (Figure 8e), and the N₂ hydrate mass fraction increased from 0.00 to 0.35 m near the injection well (Figure 8f).

4.4.3. Potential of CO₂ Storage and CH₄ Recovery

After the 306 h injection, the well was shut down and the operations forwarded to the production stage. After a 4.5-day shut-in period, depressurization was employed in the production phase. According to the operation processes detailed by [26] the CO₂-EGHRSim V1.0 simulator was used to simulate the production. Figure 9 shows the simulated and measured accumulations of CH₄, CO₂ and N₂ with the production time.

At the end of production, the calculated CH₄ quantity recovered from the reservoir was about 1.02 × 10⁴ kg (red dash-dot line in Figure 9), the quantity of produced CO₂ was about 0.11 × 10⁴ kg (green dash line Figure 9) and that of N₂ was 0.38 × 10⁴ kg (blue solid line Figure 9). In the report on the field trial by [26], the produced quantity of CH₄ was about 0.97 × 10⁴ kg (red squares in Figure 9), that of CO₂ was 0.083 × 10⁴ kg (green deltas in Figure 9), and that of N₂ was 0.39 × 10⁴ kg (blue gradients in Figure 9), respectively. Comparing the simulated results with the measured values indicated that the simulated results matched the measured values well. The relative errors were calculated as 4.60%, 26.21% and 3.22% for CH₄, CO₂ and N₂, respectively.

To quantify the recoverable capacity, a recovery efficiency of CH₄ from the hydrate-bearing reservoirs was adopted in this work, defined by referring to the recovery efficiency of oil in the CO₂-enhanced oil recovery.

$$Re=\frac{{MCH}_{4,R}}{{MCH}_{4,IAZ}}\times 100\%$$

(25)

where Re means the recovery efficiency; MCH_4,R is the recovered CH₄ from the hydrate-bearing reservoirs, in kg; MCH_4,IAZ is the initially present CH₄ in the zones which are expected to be influenced by gas production, in kg.

To evaluate the CO₂ storage potential of the hydrate-bearing reservoirs, a storage ratio was employed:

$$Sr=\frac{{MCO}_{2,injected}-{MCO}_{2,recovered}}{{MCO}_{2,injected}}$$

(26)

where Sr is the storage ratio of CO₂; MCO_2,injected means the injected CO₂ at the end of injection, in kg; and MCO_2,recovered labels the recovered CO₂ from the reservoir at the end of production, in kg.

By summarizing the hydrates initially present in the blocks affected during the field trial, the quantity of CH₄ was 1.02 × 10⁴ kg. The recovered quantity CH₄ was about 3.93 × 10⁴ kg, thus a recovery efficiency of 25.95% was obtained. In terms of CO₂, the simulated quantity of injected CO₂ was 0.25 × 10⁴ kg, and the simulated value of produced CO₂ was 0.11 × 10⁴ kg, thus the storage ratio was about 0.58.

4.4.4. Changes of Physical Characteristics due to the Recovery of CH₄

The saturation distribution of gas phase is displayed in Figure 10a. As shown in Figure 10a, the gas mainly accumulated around the production well, and the saturation of the gas phase was lower. Due to the drop of pressure around the well, the hydrate dissociated, and the released gases were pumped out, which caused the decrement in hydrate saturation near the production well. As shown in Figure 10b, near the production well, the hydrate saturation almost dropped to zero.

The composition of gas mixture was analyzed by calculating the mass fractions of CH₄ in the gas phase in Figure 11a. We found that CH₄ dominated the gas phase with a high mass fraction. Figure 11b–d illustrate the mass fractions of CH₄ hydrate, CO₂ hydrate and N₂ hydrate in the hydrate phase, respectively. Around the production well, CH₄ almost disappeared from the hydrate phase. The mass fraction of CH₄ hydrate in the hydrate phase increased along the radial distance from the production well.

5. Conclusions

In this study, a simulator, CO₂-EGHRSim V1.0, capable of simulating the replacement of CH₄ by injecting CO₂ into hydrate-bearing reservoirs was developed and verified based on experimental data. Subsequently, the simulator was applied to the Iġnik Sikumi field trial. The following conclusions can be drawn:

(1)

A thermodynamic model for determining the stable conditions of mixed gas hydrate was coupled with Tough+Hydrate V1.5, and thus a simulator, CO₂-EGHRSim V1.0, was developed in this work. The verification of the simulator suggested a high accuracy and an acceptable reliability.

(2)

The application of CO₂-EGHRSim at the Iġnik Sikumi field indicated that during the injection of a gas mixture, a hydrate mostly formed at the plume front of the gas phase; however, the CO₂–CH₄ exchange mostly occurred inside the gas plume. At the end of production, the CO₂ hydrate was only present in the reservoir near the production well, where CH₄ was mostly recovered.

(3)

The simulated CO₂ storage ratio and CH₄ recovery efficiency were 0.58 and 25.95%, respectively. The difference of simulated storage ratio between the measured value can be explained by using a thermodynamic model to predict the stable conditions for the mixed hydrate.

This study can provide guidelines for carrying out hydrate-based CO₂ marine geological storage and for analyzing the implementation feasibility of the CH₄–CO₂ exchanging method on a field scale.