AN ANALYTICAL SOLUTION AND NUMERICAL MODELING STUDY OF GAS HYDRATE SATURATION EFFECTS ON POROSITY AND
PERMEABILITY OF POROUS MEDIA
by Fangyu Gao
c
Copyright by Fangyu Gao, 2015 All Rights Reserved
A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial fulfillment of the requirements for the degree of Master of Science (Petroleum Engineering).
Golden, Colorado Date Signed: Fangyu Gao Signed: Dr. Luis Zerpa Thesis Advisor Golden, Colorado Date Signed: Dr. Erdal Ozkan Professor and Head Department of Petroleum Engineering
ABSTRACT
A gas hydrate is a type of crystallized compound formed by small gas molecules and water under high pressure and low temperature. Natural gas hydrate reservoirs exist mostly in offshore areas of outer continental margins, and some also occur in permafrost areas. Worldwide methane content in gas hydrate accumulations has an estimated volume ranging from 500 Tcf to 1.2 million Tcf. Economic values of these gas hydrate reservoirs are tremendous. A better understanding of the properties of gas hydrate-bearing reservoirs could lead to the development of novel safe and economic production methods.
There are two major deposition types of gas hydrate in porous media: pore filling and grain contact cementing-not including those related to geomechanics effect during research of hydrate-bearing reservoirs. The difference between these two distribution types is related to a nucleation condition at the beginning of hydrate cluster formation. Previous work from Verma and Pruess, 1988, shows that hydrate distribution in the pore volume and pore throat depends on the length difference for correlation between the porosity and permeability change in porous media. The previous correlation only considers the contact cementing deposition type. This thesis focuses on the study of correlation models between permeability and porosity changes during formation and dissociation of gas hydrates in porous media, called the permeability adjustment factor. A series of equations has been developed based on consideration of different parameters in the correlation, such as the power factor and critical porosity.
In previous permeability adjustment factor equations, permeability is calculated with the perme-ability equation based on a tubes-in-series model-only considering the contact cementing deposition type. This research first derives the permeability equation for the pore filling deposition type. To give a more comprehensive equation considering both types of deposition, MATLAB has been used to correlate permeability and porosity. This resulted in equations that involve a combination of pore filling and contact cementing deposition types. Finally, the TOUGH+Hydrate, T+H, numer-ical simulator was modified to include these equations to show the different results of production and saturation results for a depressurization process.
With modified permeability adjustment factors that better represent hydrate distribution in pore spaces, the T+H simulations will be more accurate. This will contribute to the future economic production of gas hydrate.
TABLE OF CONTENTS
ABSTRACT . . . iii
LIST OF FIGURES . . . vii
LIST OF TABLES . . . x NOMENCLATURE . . . xi LIST OF ABBREVIATIONS . . . xv ACKNOWLEDGMENTS . . . xvi DEDICATION . . . xvii CHAPTER 1 INTRODUCTION . . . 1 1.1 Thesis objectives . . . 2 1.2 Thesis layout . . . 2
CHAPTER 2 CHARACTERISTICS OF GAS HYDRATES . . . 4
2.1 Classification of gas hydrates based on deposition types . . . 5
2.2 Hydrate precipitation and formation processes in porous media . . . 6
2.3 Micro structure of gas hydrates . . . 7
2.4 Hydrate nucleation process in porous media . . . 9
2.4.1 Homogeneous nucleation . . . 13
2.4.2 Heterogeneous nucleation . . . 15
2.5 Relative permeability models . . . 16
2.5.1 OPM model . . . 18
2.5.2 EPM model 1 . . . 19
CHAPTER 3 METHODOLOGIES AND APPROACHES . . . 22
3.1 Addressing the problems from previous work . . . 22
3.2 Methodologies and approaches . . . 23
3.2.1 Combined model for pore filling and contact cementing hydrate distributions 25 3.2.2 Implementation of permeability adjustment factor models in a hydrate reservoir simulator . . . 30
CHAPTER 4 RESULTS AND DISCUSSIONS . . . 35
4.1 Input parameters of hydrate reservoir simulation . . . 36
4.2 Results analysis from original code . . . 37
4.3 Tecplot results in terms of pressure and temperature distribution obtained with original code . . . 40
4.4 Results of modified permeability adjustment factor equation with critical porosity and without permeability exponent . . . 42
4.5 Results of modified permeability adjustment factor equation with neither critical porosity nor permeability exponent . . . 48
4.6 Results of modified permeability adjustment factor equation with both critical porosity and permeability exponent . . . 52
4.7 Results of modified permeability adjustment factor equation without critical porosity but with permeability exponent . . . 55
4.8 Overall discussion of the results . . . 57
CHAPTER 5 CONCLUSION . . . 61
CHAPTER 6 SUGGESTIONS FOR FUTURE WORK . . . 62
REFERENCES CITED . . . 64
APPENDIX A - DERIVATION OF PERMEABILITY IN CONCENTRIC TUBES . . . 67
APPENDIX B - DERIVATION OF PERMEABILITY ADJUSTMENT FACTOR EQUATION . . . 71
APPENDIX C - T+H INPUT FILE EXAMPLE . . . 74
LIST OF FIGURES
Figure 2.1 Worldwide gas hydrates reservoir distribution . . . 4
Figure 2.2 Classification of gas hydrate based on deposition types. . . 5
Figure 2.3 Phase diagram of system with water, gas, ice and hydrate . . . 6
Figure 2.4 Dissociation process of hydrate in porous media due to thermal treatment . . . . 7
Figure 2.5 Three types of gas hydrates structures . . . 8
Figure 2.6 Relationship between water density and temperature under standard atmospheric pressure . . . 10
Figure 2.7 Crystal formation as a function of subcooling relative to the equilibrium line (AB) and the spinodal line (CD super saturation limit). . . 12
Figure 2.8 Schematic of the formation of a critical nucleum according to Classical Nucleation Theory . . . 13
Figure 2.9 Hydrate formation as a function of sub-cooling relative to the equilibrium line (AB) and the spinodal line (CD; super-saturation limit) . . . 14
Figure 2.10 Illustration of cluster formation locations (Reproduced from ). . . 16
Figure 2.11 Example of lab data plot that is usually used for relative permeabilities and saturations . . . 17
Figure 3.1 Tube in a series effect. . . 23
Figure 3.2 Example of flow between concentric tubes. . . 25
Figure 3.3 Tubes-in-series model with random sizes of pore and throat sizes. . . 26
Figure 3.4 Grain size, pore size, and pore-throat size for 27 sandstone samples . . . 27
Figure 3.5 Weibull’s distribution with shape factor equals to 1, 2, 3 and 4. . . 27
Figure 3.6 MATLAB result with Weibull’s distribution shape factor equals to 7. . . 28
Figure 4.2 Pressure and temperature results from original code. . . 38
Figure 4.3 Relative permeability curves. . . 39
Figure 4.4 Gas and water production rate of results from original code. . . 39
Figure 4.5 Pressure change at different time steps during dissociation. . . 41
Figure 4.6 Temperature change at different time steps during dissociation. . . 41
Figure 4.7 Hydrate saturation at time step 1 and 2. . . 41
Figure 4.8 Distribution of aqueous saturation in different time steps. . . 42
Figure 4.9 Gas saturation distribution as in time step 1 and 2. . . 43
Figure 4.10 Gas production profile for modified simulation with critical porosity and no permeability exponent. . . 44
Figure 4.11 Water production profile for modified simulation with critical porosity and no permeability exponent. . . 45
Figure 4.12 Gas saturation for modified simulation with critical porosity and no permeability exponent. . . 46
Figure 4.13 Water saturation for modified simulation with critical porosity and no permeability exponent. . . 47
Figure 4.14 Hydrate saturation for modified simulation with critical porosity and no permeability exponent. . . 47
Figure 4.15 Ice saturation for modified simulation with critical porosity and no permeability exponent. . . 48
Figure 4.16 Gas production profile for modified simulation with critical porosity and no permeability exponent. . . 49
Figure 4.17 Water production profile for modified simulation with critical porosity and no permeability exponent. . . 49
Figure 4.18 Gas saturation for modified simulation with critical porosity and no permeability exponent. . . 50
Figure 4.19 Water saturation for modified simulation with critical porosity and no permeability exponent. . . 50
Figure 4.20 Hydrate saturation for modified simulation with critical porosity and no permeability exponent. . . 51
Figure 4.21 Ice saturation for modified simulation with critical porosity and no
permeability exponent. . . 51 Figure 4.22 Gas production profile for modified simulation with critical porosity and
permeability exponent. . . 53 Figure 4.23 Water production profile for modified simulation with critical porosity and
permeability exponent. . . 54 Figure 4.24 Gas saturation for modified simulation with critical porosity and permeability
exponent. . . 54 Figure 4.25 Water saturation for modified simulation with critical porosity and
permeability exponent. . . 55 Figure 4.26 Hydrate saturation for modified simulation with critical porosity and
permeability exponent. . . 56 Figure 4.27 Ice saturation for modified simulation with critical porosity and permeability
exponent. . . 56 Figure 4.28 Gas production profile for modified simulation with critical porosity and
permeability exponent. . . 57 Figure 4.29 Water production profile for modified simulation with critical porosity and
permeability exponent. . . 58 Figure 4.30 Gas saturation for modified simulation with critical porosity and permeability
exponent. . . 58 Figure 4.31 Water saturation for modified simulation with critical porosity and
permeability exponent. . . 59 Figure 4.32 Hydrate saturation for modified simulation with critical porosity and
permeability exponent. . . 59 Figure 4.33 Ice saturation for modified simulation with critical porosity and permeability
LIST OF TABLES
Table 2.1 Solution properties of natural gas hydrate components at 298 K . . . 11
Table D.1 T+H input parameters for memory section . . . 76
Table D.2 T+H input parameters for rocks section . . . 77
Table D.3 T+H input parameters for hydrate section . . . 77
Table D.4 T+H input parameters for computational section . . . 78
Table D.5 T+H input parameters for computational section continued . . . 79
NOMENCLATURE
A . . . Surface area b . . . Klinkenberg b-factor accounting for gas slippage effects CR . . . Heat capacity of the dry rock FA . . . Area adjustment factor fb,j . . . Bulk phase experimental fugacity of component j feq . . . fugacity at equilibrium at temperature T Fκ . . . Darcy flux vector of component κ Fφs . . . Permeability adjustment factor FP T . . . Porosity adjustment factor fσ . . . Radiance emittance factor fv . . . Fugacity in the gas phase at temperature T f∞,j . . . Equilibrium fugacity of component j g . . . Gravitational acceleration vector Gcritical . . . Critical Gibbs free energy Gs . . . Surface excess Gibbs free energy Gv . . . Volume excess Gibbs free energy hβ . . . specific enthalpy of phase β k . . . Permeability K0 . . . Absolute permeability at large gas pressures K0 . . . Intrinsic hydration reaction constant k0 . . . Initial permeability
k00 . . . Reference permeability KR . . . Thermal conductivity of the rock Kβ . . . Thermal conductivity of phase β krβ . . . Relative permeability of phase β krφ . . . Relative permeability affected by porosity krS . . . Permeability S-factor Mκ . . . Mass accumulation term of component κ n . . . Inward unit normal vector n . . . Simple power describing dependency of porosity on permeability nw . . . Water molecule number per gas molecule PA . . . Pressure of the aqueous phase P − P∞ . . . Overpressure qκ . . . Source and sink term of component κ r . . . Radius R . . . Universal gas constant rc . . . Critical radius SA . . . Aqueous saturation S0 . . . Solubility of a bulk crystal with effectively curvature Sβ . . . Phase Saturation of β Sh . . . Hydrate saturation Sr . . . Solubility of a spherical crystal of radius r S∗ . . . Scaled saturation t . . . Time T . . . Temperature
Uβ . . . Specific internal energy of phase β UH . . . Specific enthalpy of hydrate dissociation and formation V . . . Volume vh . . . Molar volume of hydrate Vm . . . Molar volume Vn . . . Volume of subdomain n vw . . . Molar volume of water Xβκ . . . Mass fraction of component κ in phaseβ
Greek Letters
αP . . . Pore compressibility αT . . . Pore thermal expansivity β . . . Shape factor of Weibull’s distribution ∆Ea . . . Hydration activation energy ∆G . . . Change of Gibbs free energy ∆Gcritical . . . Change of critical Gibbs free energy ∆G0critical . . . ∆Gcritical when contact angle between hydrate crystal and a surface is 0 ∆gv . . . Change of Gibbs free energy per unit volume γ . . . Surface tension Γn . . . Surface area of subdomain n κ . . . component µA . . . Viscosity of the aqueous phase φ . . . Porosity φ Fraction times ∆Gcriticalin homogeneous nucleation to obtain that in heterogeneous nucleation
φc . . . Critical porosity φrr . . . porosity relative magnitute porosity φ00 . . . Reference porosity π . . . Ratio of a circle’s circumference to its diameter ρβ . . . density of phase β ρR . . . Rock density σ . . . Surface tension σ0 . . . Stefan-Boltzmann constant (5.6687 × 10−8J m−2K−4) θj . . . fractional filling of the hydrate cages on a water free bases θ . . . Contact angle between hydrate crystal and a surface
LIST OF ABBREVIATIONS
BLNL . . . Berkeley Lawrence National Lab EPM . . . Evolving Porous Medium OOP . . . Object-Oriented Programming OPM . . . Original Porous Medium T+H . . . TOUGH+Hydrate
ACKNOWLEDGMENTS
I would like to express my special appreciation and thanks to my advisor Dr. Luis Zerpa, you have been guiding and supporting me from the very beginning of my research. I would like to thank you for opening the door of the beautiful world of gas hydrates. I would also like to thank you for trusting me and financially support me, and because of that, I was able to concentrate on my research and study with all my time and passion.
I would also like to thank my committee members, Dr. Xiaolong Yin and Dr. Carolyn Koh. I would like to thank you for your help and support in knowledge unconditionally. Your valuable comments and suggestions for my proposal and research guided me into a correct path of pursuing knowledge.
I would like to thank my friend Shihao, for generously guiding me with his knowledge.
I would like to thank Dr. Mark Miller for helping me with grammar and formatting correction of my thesis.
I would like to thank Denise and Terri for helping me with all administration works. I would like to thank JK for distracting me but teaching me valuable knowledge.
I would like to thank all my friends and family for always being there and giving me supports in the hard times.
CHAPTER 1 INTRODUCTION
In the 1960s people started finding a type of “flammable ice”. That “flammable ice” is now known as gas hydrate. Gas hydrate reservoirs are widely distributed in permafrost continents, continental slopes, uplifted active and passive continental margins, oceans, and deep lakes (Collett, 2008). Although almost every country has some gas hydrate reservoirs, most of them are offshore on outer continental margins and some also occur in permafrost in the south and north poles (Potential Gas Committee, 1980). Microscopically, gas hydrate is a compound consisting of repeated hydrate molecules. Each hydrate molecule consists of water and methane molecules, and can be seen as a basic unit cell. Under standard conditions, one unit cell of gas hydrate can dissociate to a maximum of 164 unit cells of methane gas (Sloan, 1998). Based on the Potential Gas Committees report in the 1980s, methane gas hydrate reserves in permafrost regions have 14 to 3400 trillion cubic meters, and all known gas hydrate reserves has the volume of 7.6 exa (1018) cubic meters (Mosher et al., 2005). Because of the surprising number of reserves, natural gas hydrate could be considered an important potential alternative energy resource.
Currently, depletion of existing oil reservoirs is tremendous, and finding alternative energy reserves is critical to the future development of human beings (Makogon et al., 1981). Gas hydrate is not only a type of unconventional gas reservoir but also a type of potential alternative energy. Based on seismic data and other exploration methods, gas hydrate reservoirs are distributed in the deep ocean and permafrost areas. The reservoirs in ocean regions solely have an area of around 40 million square meters, almost a quarter of the ocean area (Zheng, 2011). Up to 2009, there were 116 gas hydrate reservoirs found worldwide (McGuire et al., 2009). The volume of methane gas stored in those gas hydrate reservoirs is much larger than the equivalent known conventional or shale gas reservoirs. Although gas hydrate reservoirs bring hope to the energy shortage problem in the future, a lot of technical issues need to be solved before turning to commercial production. Gaining more knowledge and a better understanding of gas hydrate reservoir behavior and properties can help in making safe and efficient economic production of gas hydrates a reality.
Estimating the production from a gas hydrate reservoir using numerical reservoir simulation is very important. With numerical simulation, many uncertainties can be solved based on the relationship between known information and unknown factors. Reliable results can be obtained based on good understanding of each parameter and an accurate phenomenon of the process. Current reservoir simulators, such as TOUGH+Hydrate, use porosity and permeability models that are not directly representing hydrate systems in reality. This thesis is focused on the study of hydrate saturation effects on porosity and permeability of porous media, during gas hydrate formation and dissociation.
1.1 Thesis objectives
The purpose of this project is to improve current comprehension of the relationship between permeability and porosity during precipitation in hydrate-bearing sediments. Following are the steps to obtain the result.
• Review literature for the characters and properties of hydrate-bearing sediments
• Analyze procedures to get current expression of relationship between permeability and poros-ity
• Formulate the limitations and disadvantages of current correlation
• Derive equations of permeability for flows between two concentric cylinders
• Develop a correlation to better describe relationship between permeability and porosity for pore filling, contact cementing and combination of both hydrate distribution types using MATLAB
• Describe production profile and saturation results of the new correlations between permeabil-ity and porospermeabil-ity with TOUGH+Hydrate simulator with modified code
1.2 Thesis layout
The layout of this thesis paper follows the order of objectives above. Below are the abstracts for each chapter.
Chapter 2 Literature review of concept and background of the research
Chapter 3 Addressing previous problems, providing methodologies and approaches as the thesis of the research as well as a solution to the previous problem
Chapter 4 Showing the results and discussion Chapter 5 Conclusion of the whole thesis
CHAPTER 2
CHARACTERISTICS OF GAS HYDRATES
Gas hydrate reservoirs are found mostly in sediments under permafrost areas and the ocean floor. Worldwide distribution of gas hydrate reservoir locations is shown in Figure 2.1. Conditions of gas hydrates formation are high pressure, low temperature, and the presence of natural gas and water. Natural conditions in these two types of areas match those of forming gas hydrate. Natural gas hydrates accumulations are always found from the drilling process or seismic data. However, in order to get the precise location and resource estimate for gas hydrates, both drilling and seismic interpretation needs to be done in the candidate area (Vanneste et al., 2001).
Figure 2.1: Worldwide gas hydrates reservoir distribution (USGS, 2012).
After identifying the location of gas hydrate reservoirs, further reservoir modeling and micro-scale simulation will help estimate a more accurate reserve. Calculation for gas hydrates reservoir is similar to that of normal conventional gas reservoirs, except for the variable solid phase saturation. As in estimation of conventional reserves, gas hydrate reserve is calculated with rock properties
volume factor (Moridis et al., 2008). However, the uncertainty of porosity and permeability makes estimating gas hydrate reserves hard.
Because solid crystal-like phase forms or dissociates with changes of pressure and temperature in gas hydrate reservoirs, permeability and porosity of the reservoir do not remain constant with time. Gas hydrates stability can be presented in a diagram as a dependence of pressure and temperature. Based on such charts we can estimate the stable region of solid gas hydrates in the reservoir under different conditions. Once we know the porosity of the gas hydrates reservoir under one specific pressure and temperature, we can use correlations to estimate permeability and reserve. However, the reality is not as simple as one single correlation equation. The truth is permeability is also dependent on pore channel geometry and deposition occurrence process of a solid in porous sediments.
2.1 Classification of gas hydrates based on deposition types
There are three classes of natural gas hydrate accumulation. Class 1 has an overlying hydrate bearing interval and a two-phase fluid interval with existence of mobile gas (Moridis, 2004). Class 2 has a hydrate-bearing overlying zone and a mobile water zone. Class 3 has only a single hydrate zone, and there is no mobile fluid zone present. Class 1, 2 and 3 of gas hydrates deposition types are shown in Figure 2.2.
Figure 2.2: Classification of gas hydrate based on deposition types.
Class 1 is the ideal target production deposition type because it has guaranteed gas production. Also, due to the coexistence of gas, water, and hydrates in the fluid interval, the hydrate zone
interface is at thermodynamic equilibrium with the two phase zone. Under this condition, a small change in temperature or pressure will dissociate the hydrate. A Class 2 hydrate reservoir requires more depressurization, thermal or inhibition treatments to release gas from hydrate, and Class 3 is the most difficult one to have gas production (Moridis et al., 2007).
Distribution of gas hydrates is related to the deposition type, and the specific mechanisms are presented in Section 2.3 and 2.4.
2.2 Hydrate precipitation and formation processes in porous media
Gas hydrates have the tendency to precipitate in the low-temperature and high-pressure zone first, and dissociate in the high temperature and low pressure zone.
Conditions of hydrate equilibrium can be described as a function of pressure and temperature, and the phase diagram is shown in Figure 2.3. In the phase diagram, I represents for ice, H is hydrate phase, Lw is liquid water, and V is gas. It is hard to see pure solid gas hydrate exist naturally, so in the phase diagram, hydrates and ice coexist in the high pressure low temperature zone.
There are three major ways to dissociate gas hydrate in order to produce gas. These three ways are depressurization, thermal stimulation, and chemical inhibition. Depressurization method uses decreasing reservoir pressure to make gas hydrate in the reservoir unstable and release methane molecules from the water cage. The thermal method is described below, and it is basically increasing temperature to dissociate hydrate. The third way is using chemical inhibitors to prevent the formation of solid hydrate. It is used mostly during the transportation of gas after production. Popular inhibitors that are used to prevent the formation of gas hydrates include alcohols and electrolytes.
Figure 2.4 explains the thermal method of dissociating gas hydrates to produce methane from the formation. When hot water is injected into the gas hydrate zone, the temperature in the near wellbore location increases because of heat transfer. Gas hydrates are no longer stable because although pressure remains the same, temperature goes up. As a result, gas hydrates near the wellbore start to dissociate. After dissociation, gas hydrate is no longer in solid form, and free gas and water are pushed into the wellbore and produced to the surface.
Figure 2.4: Dissociation process of hydrate in porous media due to thermal treatment (Lehmkoster et al., 2014).
2.3 Micro structure of gas hydrates
Gas Hydrates are a type of compound with a natural gas molecule being trapped inside a hydrogen-bonded water cage. Under low temperature and high pressure, when the compound is crystallized from water and methane gas, the resulting crystals are hydrates. Based on the formation type of gas hydrate, we can see that gas hydrate is not simply the mixture of gas and water. Indeed,
the formation of gas hydrate is a chemical reaction, and gas hydrate is a pure substance.
Each unit cell of gas hydrate crystal contains an identified number of water molecules, and the number is determined by structure type of gas hydrates. There are three main types of structure for gas hydrates: Type I, Type II, and Type H. Type I and Type II are the common types and Type H is rare to observe. Those three types are shown in Figure 2.5. In this paper, only type I will be considered since it is only related to methane and light hydrocarbons.
Figure 2.5: Three types of gas hydrates structures (Sloan, 2003).
One Type I gas hydrate unit cell has two types of cages: 512and 51262. 512has a radius of 0.395 nm, and it is called the small cage. It has a shape of pentagonal dodecahedron, which contains 12 pentagons in the cage. 512 62 cage has the radius of 0.433 nm, and it is called the large cage. It is a hexagonal truncated trapezohedron, which consists of 12 pentagons and 2 hexagons. Each Type I gas hydrate unit cell consists two 512 cages and six 512 62 cages with an average radius of 1.2 nm, and has 46 water molecules. Besides methane, Type I gas hydrate cage can also trap carbon dioxide as a common guest molecule.
A Type II unit cells has two types of cages: 512 and 512 64. 512 cage in Type II has the same shape as in Type I, but it has a smaller radius of 0.391 nm, and it is the small cage. 512 64 is the large cage in Type II. It is a hexa-decahedron, which means it contains 12 pentagons and 4 hexagons. 512 64 cage has a radius of 0.473 nm (von Stackelberg and Muller, 1954). Each Type II unit cell contains sixteen 512 cages and eight 512 64 cages. One Type II unit cell has a radius of 1.73 nm, and it contains 136 water molecules. Normally, a Type II gas hydrate cage traps oxygen and nitrogen in its common formation.
Type H unit cell has three basic shapes of cages. As the other two types, it has 512 cage with an average radius of 0.391 nm. It also contains two other types of cages: 43 56 63 and 512 68 . 43 56 63 cage is a small cage with three quadrilateral, 6 pentagons, and 3 hexagons. 512 68 cage is a very large cage with 12 pentagons and 8 hexagons. One type H unit cell has three 512 small cages, two small 43 56 63 cages, and one very large 512 68 cage. It has a radius of 1.226 nm and 34 water molecules. Type H gas hydrate always trap large molecules such as butane because of its large structure (Sloan, 1998). It is frequently observed in the Gulf of Mexico because supplies of heavy hydrocarbons produced thermogenically are common there.
2.4 Hydrate nucleation process in porous media
Nucleation of gas hydrate is the first step of hydrate formation process. What makes hydrate nucleation important and difficult to simulate is its stochastic characterization. During the nu-cleation process, tens to thousands of water and gas molecules move randomly in the fluid body (Mullin, 1993). When random motion of those molecules makes a cluster and treaches the mini-mum cluster size of hydrate molecule, the initial cluster of gas hydrate forms. This process is called nucleation, and this cluster is the beginning of solid gas hydrate. The process of nucleation is very hard to predict in lab-level experiments because of the micro-scale phenomenon and the fast motion of water and gas molecules. As a result, nucleation of gas hydrate is statistically probable but not deterministically certain (Sloan, 2003).
Because of the uncertainty of the nucleation process, a lot of hypothesis of gas hydrates are made based on our understanding of ice. Also, as a result, the TOUGH+Hydrate (T+H) simulator ignores the pore filling deposition type, since it is an “unlike-to-achieve” circumstance. However, the simulator should provide a more realistic solution instead of simply assuming one of the deposition
types doesn’t occur at all. To accomplish this task, we need to have a better understanding on the important factors that will affect the initial cluster appearance locations and hydrate nucleation types. There are three important concepts that will affect the initial cluster appearance locations: supercooling water, the solubility of gas in water, and meta-stability of ice and hydrate.
Supercooled water is a special condition in the phase behavior. It remains liquid below its freezing point. Under this condition, water is meta-stable. Water has some special characteristics compared to other types of liquid fluid. For example, the density of water is highest when the temperature is close to 4 degree Celsius. It decreases when the temperature is higher or lower than 4 degree Celsius under standard atmospheric pressure. The relationship is shown in Figure 2.6. Also, volume and entropy fluctuations of liquid water are larger when temperature drops. This feature of water helps to form an open hydrogen bonded network, and will affect the thermodynamics of water. As a result, volume increases with decreasing entropy, (Debenedetti, 2003).
Figure 2.6: Relationship between water density and temperature under standard atmospheric pressure (Stevens, 2014).
Formation of hydrate is similar to that of ice, and the supercooled water concept can be used to understand the pre-cluster stage of hydrate (Makogon et al., 1981). As shown in lab experiments, the key points of the existence of highly super-cooled water are the characters of hydrogen-bonded polyhedral. There are two main characters: concentration and location distribution (Ito et al., 1999; Speedy, 1987). These polyhedrons are inside the hydrogen-bonded molecular cages. The angles of hydrogen-bonds make the polyhedral likely to share edges and faces, thus more likely to attract each other. This is the important factor of the pre-cluster stage for the hydrate. This effect is also known as “cage effect” (Chen et al., 1997).
The solubility of gas in water is always low for non-polar gasses. Solubility and thermodynamics of different gasses are shown inTable 2.1. From this table, we can understand why reservoirs with different non-polar gas components act differently, and the reasons for clusters that can form in non-polar molecules (Sloan and Koh, 2007). The tendency of water clusters with dissolved gas to share edges and faces is the same as that for pure super-cooled water as mentioned in previous paragraph (Stillinger and Weber, 1981).
Table 2.1: Solution properties of natural gas hydrate components at 298 K
Component Solubility [105x2] −∆Hsoln [KJ/mol] −∆Ssoln[J/(Kmol)] ∆Cp [kJ/(Kmol)]
Methane 2.48 13.26 44.5 55 Ethaqne 3.10 16.99 57.0 66 Propane 2.73 21.17 71.0 70 Iso-butane 1.69 25.87 86.8 NA Nitrogen 1.19 10.46 35.1 112 Hydrogen solfide 177.9 526.35 88.4 36 Carbon dioxide 60.8 19.43 65.2 34
Normally, only line AB is shown in the diagram of crystallization, which demonstrates the relationship between concentration and temperature in Figure 2.7. This means once temperature and concentration are on line AB or on the left side of AB, the crystal will form. However, there is actually a meta-stable zone in which crystal is likely to form, but not necessarily. Furthermore, line CD is the line of super-saturation. Most of the time, if the solution is under the concentration and temperature conditions on the left side of line CD, crystallization occurs because of the force caused by supersaturation condition of the solution (Mullin and Jancic, 1979).
Figure 2.7: Crystal formation as a function of subcooling relative to the equilibrium line (AB) and the spinodal line (CD super saturation limit)(Sloan and Koh, 2007).
Figure 2.8 shows the initial growth of a cluster. Nucleation of gas hydrates is the first step of the gas hydrate formation. This process needs water molecules to combine with gas molecules to form a sufficient size of the cluster to break the surface energy and form a new contact interface. In order to obtain enough energy and form the initial small cluster, a super-saturation condition is necessary (Clennell, 1991). Kelvin equation shown below in Equation 2.1 describes the relationship between super-saturation and the initial size of clusters that can form gas hydrates.
ln(S r S0) = 2γ r ¯ Vm RT (2.1)
In this equation, Sr and S0 represents the solubility of a spherical crystal of radius r and a bulk crystal with effectively-curvature respectively. While γ is the specific surface energy, r is the nucleating crystal radius, and Vm is molar volume. Supersaturation is a result of decreasing temperature, increasing pressure, or increasing the concentration of methane (Clennell, 1991). In a regular phase diagram of a system with water, gas, ice and hydrate, when it rehearses conditions of super-saturation, it is already in hydrates zone.
Figure 2.8: Schematic of the formation of a critical nucleum according to Classical Nucleation Theory (Hester et al., 2006).
Figure 2.9 is the super-saturation condition for hydrate formation, which has the same principle as crystallization of the salt in water solution. As a conclusion for Figure 2.9, the right side of line AB is the region in which no hydrate can form; the left side of line AB and right of line CD is the metastable zone where hydrate may form; left of line CD is the unstable zone where hydrate is very likely to form. However, the problem is to make the probable phenomenon of crystallization predictable and to quantify precipitation of hydrate in this area.
Although the Kelvin equation can describe conditions for the first sufficient cluster of initial hydrate nucleation, it is almost impossible to simulate the random movement before and after cluster formation. Formation criterias such as temperature and pressure are not simply a line, but a region shown in Figure 2.9. The cluster can form anywhere between the equilibrium line and the spinodal or supercooling line. In other words, it is unrealistic to include the Kelvin equation in a simulator, and finding another methodology is necessary to simulate the beginning of gas hydrate formation.
2.4.1 Homogeneous nucleation
Homogeneous nucleation is not considered in T+H and some other simulators because it has rarely happened in lab experiments, thus it is hard to observe. It is not only difficult to generate in the lab but also rare to find in real hydrate reservoirs. The reason homogeneous nucleation is hard to reach is because of the force it needs to form an initial cluster in formation fluid with dissolved
Figure 2.9: Hydrate formation as a function of sub-cooling relative to the equilibrium line (AB) and the spinodal line (CD; super-saturation limit) (Sloan and Koh, 2007).
gas. The initial energy it needs to break to form nucleation inside the fluid is much larger than the energy it needs to break for heterogeneous nucleation, which means forming hydrate on a surface or with impure fluid. The other reason for the rare appearance could then be the impurity of the formation fluid. There is always dirt, sand or other suspended substances.
Homogeneous nucleation, or in another form, pore filling hydrate formation can be formed when the formation fluid is ultra-pure. The cluster can form inside the fluid. Additionally, as Mullin (Mullin, 1993) indicates, there are 106 particles per cubic centimeter in lab condition experiments; even with carefully filtered water, there are still 103 particles per cubic centimeter. Those particles are very likely to trigger a heterogeneous nucleation.
Homogeneous nucleation requires high-level purity of the water solution with dissolved gas. The nucleation can only occur when the random motion of molecules reaches the critical size of the cluster which is energetically favorable. The energy is also called Gibbs free energy, and it is the summation of the surface excess Gibbs free energy and the volume excess Gibbs free energy. The
equation can be expressed as (Sloan and Koh, 2007): ∆G = ∆Gs+ Gv = 4πr2σ + 4 3πr 3∆g v (2.2)
Critical cluster size can be expressed as:
rc= − 2σ ∆gv
(2.3) Critical Gibbs free energy to form cluster is:
∆Gcritical= 4πσr2c
3 (2.4)
In addition, Equation 2.5 and 2.6 are another form of critical cluster size and critical Gibbs free energy that expressed with Gibbs free energy per unit volume of hydrate formed (Englezos and Bishnoi, 1988; Englezos et al., 1987).
rc= − 2σ ∆gv (2.5) −∆gv = RT vh hX2 1 θjln fb,j f∞,j +nwvw(P − P∞) RT i (2.6) 2.4.2 Heterogeneous nucleation
Unlike homogeneous nucleation, heterogeneous nucleation is easier to achieve. Heterogeneous nucleation always requires an impure solution because it needs a foreign substance to form. This foreign substance can be dust, dirt, sand, surface of the porous media wall, or surface of a pipe.
When heterogeneous nucleation forms with a surface, there are some requirements for the contact angle (φ). Then, a new Gibbs free energy level needs to be achieved. The new equation of critical Gibbs free energy as a function of contact angle is:
∆G0critical = φ∆Gcritical (2.7)
In Equation 2.7, φ is a fraction factor that can be expressed as: φ = [(2 + cosθ)(1 − cosθ)
2]
4 (2.8)
Figure 2.10 is the illustration of cluster formation locations. Part (a) of Figure 2.10 is the spherical cluster formed inside the solution, and it is for homogeneous nucleation. Part (b) of Figure 2.10 is the cap-shaped cluster formation for heterogeneous nucleation. Part (c) of Figure 2.10
is the lens-shaped cluster formation for heterogeneous nucleation. Among those cases, part (a) is a special case of non-wetting case and could also be heterogeneous if the solution is not pure. The more hydrate wetting condition is (which means small wetting angle), the more likely it will be heterogeneous nucleation. This is because thermodynamically, the hetero cap-shaped cluster forms on the hydrate-wetting surface as heterogeneous cluster instead of a homogeneous one (Kashchiev and Firoozabadi, 2002)
Figure 2.10: Illustration of cluster formation locations (Reproduced from (Kashchiev and Firooz-abadi, 2002)).
The heterogeneous nucleation is also possible to form in between two fluid contacts, which shows the free gas in the reservoir. This is because that the layer between two solutions is more likely to be supersaturated, and has a stronger tendency to unstable zone to form hydrates nucleation. As a result, initial formation and growth of hydrate is more favorable to start from this contact layer instead of the bulk. This explains why heterogeneous nucleation is easier to form than homogeneous nucleation.
2.5 Relative permeability models
perme-occupied by that phase. The relationship can be expressed as:
krβ = krβ(Sβ) (2.9)
Figure 2.11 is an example of lab data plot that is usually used for relative permeabilities and saturations.
Figure 2.11: Example of lab data plot that is usually used for relative permeabilities and saturations (Moridis et al., 2008).
When hydrate forms in the formation, there is new solid precipitation in the pores. As a result, the relative permeability changes during this process. However, the new relative permeability is hard to predict because the hydrate precipitations are not easy to locate. Though the method based on current understanding and technology is not perfect, it is still valuable to conduct the research by creating two relative permeability models for hydrate reservoir systems.
These two models are used for the description of wettability processes, which includes relative permeability and capillary pressure in reservoir systems with the present of ice and hydrate solids. The first model is called the OPM model, which refers to original porous medium. In this model,
original porosity and intrinsic permeability remain unchanged with a new formation of hydrate in the reservoir, and fluid flow is controlled by phase saturations. The second model is called the EPM model which refers to evolving porous medium model. In this model, intrinsic permeability and porosity are changed with time because this model treats hydrate as solid which performs as part of the porous medium. The pore space is assumed to be occupied only by water and gas. There are two variations of the EPM model, which are EPM 1 and EPM 2, and will be introduced and discussed later in Sections 2.5.1, 2.5.2 and 2.5.3.
2.5.1 OPM model
OPM model is “original porous medium” model. In this model, the relative permeability of water is only a function of water saturation. Because of this feature, the OPM model is simple and does not need other information to calculate the relative permeability. When new solid ice or gas hydrates form in the system, relative permeability of gas or water still remains a function of liquid saturation. In other words, no matter how gas, ice and hydrate divide the rest of the pore, water saturation is the only factor that will affect change of relative permeability of water.
For example, in the OPM model, when solid forms in the system, the saturation of solid phases such as ice saturation and hydrate saturation increases. As a result, gas saturation and water sat-uration decreases. However, the calculation of the modified water satsat-uration is simply considering the new water saturation in the system, and ignoring the effect of new formed solid in the system. Moreover, this OPM model is assuming hydrate or ice precipitate in pores with either only water or only gas present.
In OPM model, permeability adjustment factor (FφS) can be calculated as a product of relative permeability affected by porosity and permeability S-factor which is affected by newly formed solid on intrinsic permeability (Moridis et al., 2008). The equation of the permeability adjustment factor is,
FφS = krφkrS (2.10)
Because in the OPM model, solid formed in the system does not have an influence on relative permeability, the permeability S-factor always equals one. Relative permeability affected by poros-ity equals to 1 if there is no effect from porosporos-ity on permeabilporos-ity or the effect is negligible. The
effect of porosity on relative permeability can be expressed as Equation 2.11 if there is an effect.
exp[γ(FP T − 1)] (2.11)
In this expression, gamma is an empirical parameter (Rutqvist and Tsang, 2002), FP T is a porosity adjustment factor as a function of pressure and temperature. It can be calculated as shown in Equation 2.12, in which αp and αt represent for pore compressibility and thermal expansivity. φrr is φ relative magnitude, φ00 is reference porosity.
FP T = φ φrrφ00 = φ φ0 = exp[αP∆P + αT∆T ] = 1 + αP∆P + αT∆T (2.12) Capillary pressure in OPM model is expressed as:
Pcap= sqrt( φrr krr FP T krφ ) ∗ Pcap,00(S∗) (2.13)
In this equation, the assumed relationship between φ and φ00 is as shown in Equation 2.14 which is derived from Equation 2.13, and relationship between k and k00 is shown in equation2.15 which is derived from below equation:
φ φ00 = φrrFP T (2.14) k00 k = 1 krrkrφ (2.15) kβ = kkrβwherek = k0FφS = k00krrFφS (2.16) S∗ is scaled saturation from Equation 2.17, and it is shown below:
S∗= SA SA+ SG
(2.17) Because of the limitations in the OPM model, it becomes the least favorable model to calculate relative permeability when solid forms. Therefore, two EPM models are developed to solve the limitations of OPM model. These two models are EPM 1 and EPM 2, and will be introduced in the following sections.
2.5.2 EPM model 1
EPM models stand for “evolving porous medium” models. Both EPM models solve the limi-tation of OPM model, which allow relative permeability to be updated with the formation of new solid. EPM models account influence from changes of solid saturation for liquid phase relative
permeability.
Permeability cannot be calculated simply as shown in Equation 2.18 because the relative per-meability of water and gas depend on water and gas saturation, and cannot be summed together.
k = k0(krA+ krG) (2.18)
Solid saturation is not the only factor in this case, and saturation of water and gas needs to be considered. To solve this problem, each pore in the system is regarded to be filled with only one phase, which means each pore is filled by either water only, or gas only for the fluid phase. This makes the equation of permeability S-factor become Equation 2.19 or 2.20, which is a function of solid saturation, and for the overall permeability S-factor,
krS = krA(SA= 1 − SS) (2.19)
krS = krG(SG= 1 − SS) (2.20)
if assume gas saturation and water saturation equals to each other, then the estimation becomes: krS = krA(SA= 1 − SS) + krG(SG= 1 − SS) (2.21) Relative permeability of water and gas, porosity adjustment factor, and porosity effect on relative permeability has the same calculation as described in OPM model. Capillary pressure for EPM model can be calculated with:
Pcap= sqrt( φrr krr FP T(1 − SS) krφkrS )Pcap,00(S∗) (2.22)
EPM 1 model still has some limitations because the averaging process of the permeability S-factor is not accurate, and it is only taking half and half distribution to gas and water saturation. EPM 2 model in the following section gives an alternative solution for this limitation.
2.5.3 EPM model 2
EPM 2 model is almost the same as EPM 1 model. The difference between them is that EPM 2 model has a different equation to calculate permeability S-factor. In this model, the concept of critical porosity is introduced, and the permeability S-factor equation is:
krS=
hφ0(1 − SS) − φc φ0− φc
in
Critical porosity exists when solid forms, some of the porosity gets trapped, and cannot be reached by the fluid. φ0 is the original porosity of the formation. n is an empirical term that ranges between 2 to 3 most of the time to describe the dependency of permeability on porosity. All other terms are exactly same as in EPM 1 model.
In the T+H simulator, EPM 2 model is always used. In simulation part of this thesis, EPM 2 model is the only one can be used, since the permeability adjustment factor equation derived is only related to EPM 2 model.
Although the EPM models are evolutions in the relative permeability calculation in the gas hydrate area compared to OPM model, a better model can be developed in the future, if there is a more accurate calculation of relative permeability available. The future relative permeability model for hydrate and ice should be based on theoretical assumptions, lab experiment results and field case data, and should be capable of describing the real situation of hydrate-bearing reservoirs.
CHAPTER 3
METHODOLOGIES AND APPROACHES
Solid appears in the porous media when gas hydrate forms. The new solid has the ability to block fluid that will pass through the pore by decreasing pore volume. As we know, porosity can be described as the fraction of pore volume over bulk volume. In a certain size of reservoir, when the total volume is constant, the formation of solid gas hydrates decreases the pore volume while dissociation of gas hydrates increases the porosity of the target formation.
Permeability is “capability of a porous rock or sediment to permit the flow of fluids through its pore spaces”. Normally, permeability does not necessarily have a correlation with porosity in a reservoir. For example, sandstone and shale with same porosity and grain shape will have different permeability. The reason is because the same shape of grains and packing type will give same porosity no matter how big the grains are, and the condition is assuming there are no cementing or other factors that will affect porosity. However, fluid tends to flow faster in larger pore space, and that causes a difference in permeability even the porosity is the same (Nurhono et al., 2012).
However, for a certain geometry of pore channels, permeability can be correlated with porosity. When gas hydrate forms, there could be cementation that blocks the path-way of fluid flow, and at the same time reduces the pore volume. As a result, there will be a correlation between changes in porosity and permeability in porous media during formation or dissociation of gas hydrates. 3.1 Addressing the problems from previous work
There are two equations describing the correlation between permeability and porosity. One is the simplified equation which is:
k k0 = Fφs = ( φ φ0 )n (3.1)
and the other one is the Verma and Pruess equation with critical porosity shown as: k k0 = Fφs = ( φ − φc φ0− φc )n (3.2)
Fφsis the permeability adjustment factor used to express permeability change as a dependence of porosity change during solid precipitation. In this equation, k is the permeability of the formation, n is the simple power in the range of 2 to 3. The number between 2 to 3 shows different dependency levels of permeability on porosity. If the solid gas hydrate is formed in the pore throat or in contact areas between pores, permeability will change more significantly with porosity change. In this case n will be closer to 3 due to the significance of rapid drop of permeability. The phenomenon of permeability rapid drop due to the formation of gas hydrate in pore throat is observed in both lab and field studies and production data (Moridis et al., 2008).
Figure 3.1: Tube in a series effect.
The equation with critical porosity, Equation 3.2, is applicable when pore throat or contact areas between grains are blocked by gas hydrate or ice solid. When solid forms in porous media in cementing mode, most of the pores are empty and contacts between grains are cemented by gas hydrates or ice. The behavior of pores can be described as tubes in series (Verma and Pruess, 1988) as shown in Figure 3.1.
3.2 Methodologies and approaches
The purpose of this thesis is to modify the permeability adjustment factor equation to obtain expressions that better represent the phenomena of hydrate formation or dissociation in porous media and its effect on porosity and permeability. The permeability adjustment factor is the ratio
of current permeability to initial permeability, and it is given by: k k0 = Fφs = ( φ − φc φ0− φc )n (3.3)
This equation was developed by Verma and Pruess (Verma and Pruess, 1988) based on solid pre-cipitation. In this equation, k is current permeability with new hydrate deposition, k0 is original permeability with no newly formed hydrate, and it is also before solid saturation changes; Fφs is permeability adjustment factor; φ is current porosity, φc is critical porosity and φ0 is original porosity before solid saturation changes. This equation is derived based on a tubes-in-series model and considers the permeability adjustment factor as a function of current porosity, critical porosity, initial porosity, and the empirical factor n.
Despite empirical factor n, the equation expresses that the ratio of current permeability to original permeability is the same as the difference between current effective porosity and original effective porosity. In this equation, effective porosity is subtracting critical, or “unreachable”, porosity from current or original porosity.
The problem of Verma and Pruess equation is that it only includes the case of contact cementing precipitation and does not include pore filling precipitation. As mentioned in Section 2.4, pore filling precipitation of hydrate is a result of homogeneous nucleation. Although it is rare to happen in the lab and also hard to observe in nature, it still needs to be considered and included in the equation and in the simulator to give a more accurate results.
Thus, permeability equation for pore filling precipitation is obtained from the Hagen Poiseuille equation for flow between concentric tubes, and is given by,
k = R41− R4 2− (R2 1−R22)2 lnR1 R2 8R2 1 (3.4) as permeability of flow between concentric tubes which is shown in Figure 3.2. The full procedure of derivation of this equation is presented in Appendix A.
After that, the equation can be verified using an inner tube radius of zero as the condition. This reduces to a case when there is no solid inside the pore, similar to the contact cementing case. The result shows that when boundary condition applies, the result of permeability equation for pore filling precipitation is the same as that for contact cementing precipitation, which proves that the
Figure 3.2: Example of flow between concentric tubes.
After the derivation of permeability equation, a new equation for the system combining both pore filling and contact cementing precipitation models can be derived. In order to find the relation-ship for the Combined equation, MATLAB is used to generate a system with parallel tubes-in-series system to fit the result, and the procedures are shown in following section.
3.2.1 Combined model for pore filling and contact cementing hydrate distributions MATLAB was used to generate a system consisting of a certain number of random tubes connected in series and in parallel. Then, random sizes of hydrate are specified either distributed on the surface of the pore medium, or inside the pore, or both. The tubes-in-series model with random sizes of pores and throats is shown in Figure 3.3. Tubes are connected to each other for every row, and are parallel and separated from each row. Then, the permeability for each pore can be calculated using the corresponding permeability adjustment factor equation, either pore filling or grain cementing. If there is solid inside the pore, Equation 3.4 should be used. If there is no solid inside the pore, use Equation A.29 . Next calculate average permeability with harmonic averaging method which is shown as:
Harmonic mean = 1 n a1 + 1 a2 + 1 a3 + ... + 1 an (3.5)
and it will give permeability for each row. Then, sum up the permeability from each row, and that will be the overall average permeability for this system.
Figure 3.3: Tubes-in-series model with random sizes of pore and throat sizes.
The MATLAB code includes 50 tubes in a row for 50 rows. Average throat radius is 30 × 10−7 meters, and average pore radius is 76 × 10−7meters. Those numbers are from (Nelson, 2009) which are the average values for normal sandstone formations. The figure of normal pore and throat sizes for sandstone is shown in Figure 3.4.
The random number generating distribution is Weibulls distribution, which is a distribution type widely used to generate formation pore sizes. Weibull’s distribution only gives random number above 0, and the distributions with different shape factors are shown in Figure 3.5.
The shape factor used to give a normal result is 7, and the result for shape factor 7 is shown in Figure 3.6. Shape factor 400 is also used to show an extreme result, and the result is in Figure 3.7 which is similar to Figure 3.6. In both figures, hollow dots represent permeability values calculated using exact radius for pore volume and throat. Stars represent permeability values calculated
Figure 3.4: Grain size, pore size, and pore-throat size for 27 sandstone samples (Wardlaw and Cassan, 1979).
permeability adjustment factor results of the system with only contact cementing type. Blue color shows the results for the system with only pore filling deposition type. Red color represents the results for the system with half pore filling deposition type and half contact cementing type. By comparing the plot with different Weibull’s distribution shape factor, we can see that the difference between the two cases is how concentrated the values are on the real value curves, when the shape factor is increased to 400. With shape factor equals to 7, the scenario is closer to reality, and the data are generated with the random number jumping around the real values which are the curves made from circles. While the results given from shape factor equals to 400, the random numbers are generated exactly the same as real values, and thus curves generated from random numbers which are the stars are the same as the ones generated with real values which are made of hollow circles. This shows that the random numbers are reasonable, and the correlation can be verified from the MATLAB result.
Figure 3.6: MATLAB result with Weibull’s distribution shape factor equals to 7.
After several approaches, the combined model for pore filling and contact cementing hydrate distributions is obtained. Equation 3.6 represents the average permeability as function of solid phase saturation, and Equation 3.7 represents the average permeability as function of porosity.
k k0 = (1 − Sh) × Xcc+ 1 − Sh2− (1 − Sh) 2 ln(√1 ) × Xpf (3.6)
Figure 3.7: MATLAB result with Weibull’s distribution shape factor equals to 400. k k0 = φ φ0 × Xcc− φ φ0 φ φ0 + 2 + φ φ0 ln(q 1 (1− φ φ0) ) × Xpf (3.7)
In Equation 3.6 and Equation 3.7, Xcc represents for fraction of contact cementing deposition type and Xpf represents for fraction of pore filling deposition type in solid phase.
In order to give more possibilities, four models are proposed based on Equation 3.6. Other equa-tions considering either critical porosity or permeability exponents or both are shown in Equation 3.8, 3.9 and 3.10. The modified equation considering only critical porosity effect is:
k k0 = (1 − Sh− φc) × Xcc+ 1 − Sh2− φc− (1 − Sh)2 ln(√1 (Sh) ) × Xpf (3.8)
The equation considering only permeability exponent effect is: k k0 = (1 − Sh)n× Xcc+ 1 − Sh2− (1 − Sh) 2 ln(√1 (Sh) ) n × Xpf (3.9)
The equation considering both critical porosity effect and permeability exponent is: k k0 = (1 − Sh− φc)n× Xcc+ 1 − Sh2− φc− (1 − Sh)2 ln(√1 (Sh) ) n × Xpf (3.10)
3.2.2 Implementation of permeability adjustment factor models in a hydrate reservoir simulator
In order to initialize the simulation, principles of mass and energy balance equations should be studied, as well as features of the operation system and the core code of simulator being used. The simulator used in this study is TOUGH+Hydrate version 1.0 (T+H) (Moridis et al., 2008). It is a simulator developed by Berkeley Lawrence National Lab (BLNL) based on Tough 2 concepts, underlying physics, and governing equations. T+H simulator is used for simulation of hydrate-bearing reservoirs and has the ability to calculate multi-component and multi-fluid and heat flow. T+H simulator uses dynamic memory allocation which minimized storage requirements. It also uses object-oriented programming (OOP). It is using new data structures which describe objects with which the code is based, the basic object is defined through derived data types, and properties and processes are described in modules and sub-modules.
The T+H simulator has some underlying physics assumptions. For example, Darcy’s law is used, and mechanical dispersion is negligible. It also assumes hydrate has same compressibility and thermal expansivity as ice. It neglects change and movement of the geologic medium due to temperature and pressure change. Also, dissolved salts do not precipitate when water freezes. Concentration of inhibitors does not affect thermodynamic properties of the water phase fluid. The inhibitor is stable in the temperature-pressure range of the study.
T+H simulator uses the mass balance equation for all grid blocks with integral finite difference method (Pruess et al., 1999), and the equation is
d dt Z Vn MκdV = Z ΓnFκ· d ˜A + Z Vn qκdV (3.11)
In this equation, V and Vn are volume and volume of the subdomain. Mκ is mass accumulation term of component κ. A is surface area, and Γn is the surface area of subdomain n. Fκ is Darcy flux vector of component κ. n is the inward unit normal vector. qκ is the source and sink term of component κ, and t is time.
In the simulator, components κ include hydrate, water, methane, the water-soluble inhibitor which including salt or organic substance, and heat. Phases β include solid hydrate with compo-nents of methane, water for equilibrium and hydrate for the kinetic, aqueous phase with compocompo-nents
phase with water component.
Mass accumulation terms Mκ have two expressions. One is for equilibrium condition, and the other one is for the kinetic condition. The expression for equilibrium condition is:
Mκ = X
β≡A,G,I
φSβρβXβκ· κ ≡ w, m, i (3.12) The expression for kinetic condition of mass accumulation terms is:
Mκ= X
β≡A,G,H,I
φSβρβXβκ· κ ≡ w, m, h, i (3.13) In both equations 3.12 and 3.13, phi represents for porosity, ρβ is density of phase β, Sβis saturation of phase β, and Xβκ is mass fraction of component κ.
When mass accumulation is under kinetic conditions, Equation 3.2.2 is used to describe hydrate mass component and phase (Kim et al., 1987):
QH = ∂M ∂t = −K0exp ∆Ea RT FAA(feq− fv) (3.14)
In this equation, K0is the intrinsic hydration reaction constant, ∆Eais hydration activation energy, R is the universal gas constant which is 8.314 J mol−1K−1, T is temperature, FAis area adjustment factor, A is surface area participating in the reasion, feq and fv are fugacity at equilibrium at temperature T and in the gas phase at temperature T respectively.
For heat accumulation terms, the kinetic equation includes effect from rock matrix and all phases, which is:
Mθ = (1 − φ)ρRCrT + X β≡A,G,H,I
φSβρβUβ+ Qdiss (3.15) In equatino 3.15, ρR is rock density, CR is heat capacity of the dry rock, Uβ is specific internal energy of phase β, and ∆UH is the specific enthalpy of hydrate dissociation or formation. In equation3.15, Qdiss can be expressed as:
Qdiss= ∆(φρHSH∆H0) (3.16)
for the equilibrium dissociation, and it is 0 for kinetic dissociation.
Mass fluxes of water, methane, and inhibitor which include contributions from aqueous and gaseous phases can be shown as:
Fκ= X κ≡A,G
In this expression of fluxes, components κ can be water, methane, and inhibitor. Contributions of hydrate and ice to fluxes are zero. Flux from aqueous phase can be expressed as Darcy’s law:
FA= −k krAρA
µA
(5PA− ρAg) (3.18)
In this equation, k and krA are rock intrinsic permeability and relative permeability of aqueous phase respectively. µAis the viscosity of the aqueous phase. PA is the pressure of aqueous phase, and g is gravitational acceleration vector.
Mass flux of gaseous phase incorporates advection and diffusion contributions can be shown as: FGκ = −k0(1 + b PG )krGρG µG XGκ(5PG− ρGg) + JGκ (3.19) Where, k0 is absolute permeability at large gas pressures, b is Klinkenberg b-factor for gas slippage effects. krG is relative permeability of gaseous phase, and µG is the viscosity of gaseous phase. JGκ is the term used to describe the diffusive mass flux of component κ in the gaseous phase.
Flux of inhibitor dissolved in water is:
FAi = XAiFA+ JWi (3.20)
In this equation, JWi can be expressed as: JWi = −φSW φ13S 7 3 A D0ρA5 XAi = −φSW(τA)D0ρA5 XAi (3.21) In this expression, the component κ includes water and methane. D0 represents for the molecular diffusion coefficient of inhibitor in water, and τA represents for the tortuosity of aqueous phase.
Heat flux equation includes conduction, advection and radiative heat transfer, and the equation is:
Fθ = −[(1 − φ)KR+ φ(SHKH+ SIKI+ SAKA+ SGKG)] 5 T + fσσ + 0 5 T4+ X β≡A,G
hβFβ (3.22)
Where in the equation, KR and Kβ are the thermal conductivity of the rock and different phases. hβ is specific enthalpy of the phase β, fσ is radiance emittance factor, and σ0 is Stefan-Boltzmann constant which is 5.6687 × 10−8J m−2K−4.
Mass production rate can be determined in sinks term, and injection rate can be determined in sources term: ˆ qκ = X κ≡A,G Xβκqβ (3.23)
In this expression, qβ is production rate of phase β. The occurrence of inhibitor injection can be a rate as individual which is ˆqi or as a fraction which is ˆqiXAiqˆA. In the kinetic model, the sink and source terms corresponding to hydrate dissociation and release of methane and water must be accounted for. The expression for methane production rate is:
ˆ
qm = Qm= −W m
WcQH (3.24)
Although there are three methods to use for gas hydrate dissociation simulation, only depres-surization will be studied for this research. Understanding of the T+H started from getting familiar with the user manual and sample examples. A three-dimensional kinetic dissociation gas hydrate reservoir with multiple time steps is simulated.
The target modified code is located in T M edia P roperties.f 90 file, and the section that con-tains the permeability adjustment factor equation is Adjustment-Effect of solid phases in EPM models. After modification, the files can be compiled and run successfully with Microsoft Visual Studio.
The original code is shown below:
CASE( 2 , 4 ) ! . . . . EPM Model #2 !
S f r e e = 1 . 0 d0 − S s o l i d − media ( nmat )%C r i t S a t !
IF ( S f r e e > 0 . 0 d0 ) THEN
P e r m e a b i l i t y F a c t o r = ( S f r e e / ( 1 . 0 d0 − media ( nmat )%C r i t S a t ) ) ∗∗ media ( nmat )% PermExpon
ELSE
P e r m e a b i l i t y F a c t o r = 0 . 0 d0 END IF
!
One of the modified code is shown as below. This example is with the modified equation with consideration of critical porosity and permeability exponent, and the ratio between contact cementation and pore filling deposition type is 1.0.
CASE( 2 , 4 ) ! . . . . EPM Model #2 !
S f r e e = 1 . 0 d0 − S s o l i d −media ( nmat )%C r i t S a t ! !
IF ( S f r e e > 0 . 0 d0 ) THEN
P e r m e a b i l i t y F a c t o r = ( S f r e e ∗∗ media ( nmat )%PermExpon ) ∗ 1 . 0 d0 + ( ( 1 . 0 d0 − ( 1 . 0 d0− S f r e e ) ∗∗2− S f r e e ∗∗2/ l o g ( 1 . 0 d0 / ( 1 . 0 d0−S f r e e ) ) ) ∗∗ media ( nmat )%PermExpon ) ∗ ( 1 . 0 d0 −1.0 d0 ) ! f a n g y u ELSE P e r m e a b i l i t y F a c t o r = 0 . 0 d0 END IF !
As mentioned previously, because of lacking of field and lab data, this research remains in theoretical level, and further suggestions for future research area and verification steps are in the Future Work chapter. T+H simulation helps provides more possibilities and a new direction for a better and a more accurate future simulator.
CHAPTER 4
RESULTS AND DISCUSSIONS
The results for production rate of water and gas as well as saturation for gas, water, hydrate, and ice are presented in this section.
Those results are simulated with TOUGH+Hydrate 1.0 original and modified versions. The permeability adjustment factor equation is modified into four models. Those four models are based on the correlation between porosity change and permeability change considering pore filling and contact cementing hydrate distributions. All models are based on the simple correlation of permeability change as function of porosity change with hydrate saturation:
k k0 = (1 − Sh)n× Xcc+ 1 − Sh2− (1 − Sh) 2 ln(√1 (Sh) ) n × Xpf (4.1)
which was introduced in Section 3.2.1 and its derivation is presented in Appendix A and Ap-pendix B. This equation is then modified to obtain four models based on the previous permeability adjustment factor equation derived by Verma and Pruess, 1988, which is:
k k0 = Fφs = ( φ − φc φ0− φc )n (4.2)
Those four models are: the equation of permeability adjustment factor equation with critical porosity and without permeability exponent:
k k0 = (1 − Sh− φc) × Xcc+ 1 − Sh2− φc− (1 − Sh)2 ln(√1 (Sh) ) × Xpf (4.3)
the equation with both effect from critical porosity and permeability exponent: k k0 = (1 − Sh− φc)n× Xcc+ 1 − Sh2− φc− (1 − Sh)2 ln(√1 (Sh) ) n × Xpf (4.4)
the equation with neither critical porosity or permeability exponent: k k0 = (1 − Sh) × Xcc+ 1 − Sh2− (1 − Sh) 2 ln(√1 (Sh) ) × Xpf (4.5)
and the equation without critical porosity but with permeability exponent: k k0 = (1 − Sh)n× Xcc+ 1 − Sh2− (1 − Sh) 2 ln(√1 (Sh) ) n × Xpf (4.6)
