Gas Saturation Resulting from Methane Hydrate Dissociation in a

Article pubs.acs.org/JPCC

Gas Saturation Resulting from Methane Hydrate Dissociation in a Porous Medium: Comparison between Analytical and Pore-Network Results Ioannis N. Tsimpanogiannis* Environmental Research Laboratory, National Center for Scientific Research “Demokritos”, Aghia Paraskevi, Athens 15310, Greece

Peter C. Lichtner Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States S Supporting Information *

ABSTRACT: We develop predictive tools for methane gas saturation that results from hydrate dissociation in porous media at different length scales, ranging from the single-pore scale up to the pore-network scale. Initially, we examine the case of a single spherical hydrate grain dissociating inside a bulk continuum (i.e., without the constraints from the solid surfaces that are present inside a porous medium). The growth of the resulting gas bubble is limited only by the liquid pressure of the domain, the capillary pressure of the growing bubble, and the gas solubility in the liquid phase. This case corresponds to an upper limit for the gas saturation. Next, we consider the case of a hydrate grain located inside a single pore body as well as the case of multiple hydrate grains that are distributed randomly inside the pore network. To this purpose we consider a simple porous domain that is represented by a pore network with all pores/throats being of the same size. When the hydrate phase is confined inside a porous domain, the growth of the resulting gas bubble is controlled, mostly, by the capillary thresholds of the interconnections (i.e., the pore throats) between the different pore bodies. For all cases we develop analytical solutions for the ratio of the gas to the hydrate saturation, Sg/SH, and compare the solutions with results obtained from pore network simulations [Phys. Rev. E 2006, 74, 056303] where the pore bodies/throats follow appropriate size distributions. Very good agreement is observed between the different approaches considered in this study. 12 pentagons and four hexagons (51264). Each unit cell of sII hydrate is made up of 16 small and eight large cavities.1 The unit cell of the sH hydrate (hexagonal P6mmm space group) consists of 34 water molecules that form three types of cavities: the small, a pentagonal dodecahedron (512), the medium that is formed by three tetragons, six pentagons, and three hexagons (435663), and the large that is formed by 12 pentagons and eight hexagons (51268). Each unit cell of sH hydrate is made up of three small, two medium, and one large cavity.1 In recent years a significant effort has been put into the theoretical and computational study of clathrate hydrates in order to obtain a better understanding of hydrate behavior. To this purpose, approaches such as molecular dynamics,2,3 Monte Carlo, 4,5 DFT, 6,7 ab initio,8−10 and quantum dynamic calculations11 have been utilized. The capability of clathrate hydrates to store large amounts of gases has attracted significant attention from the scientific and

1. INTRODUCTION Under appropriate conditions of high pressures (e.g., such those found in suboceanic sediments) or low temperatures (e.g., such those found in permafrost regions) water molecules can self-assemble to form three-dimensional cavities that are further stabilized by the inclusion of guest molecules. The resulting nonstoichiometric, crystalline material is known as clathrate hydrate. The hydrogen-bonded water molecules that form the hydrate lattice interact with the guest molecules with van der Waals-type forces. Based on their crystal structure, three types of hydrate structures are the most common. These structures have been identified as structures sI, sII, and sH. The unit cell of the sI hydrate (cubic Pm3n space group) consists of 46 water molecules that form two types of cavities: the small, a pentagonal dodecahedron (512), and the large that is formed by 12 pentagons and two hexagons (51262). Each unit cell of sI hydrate is made up of two small and six large cavities.1 The unit cell of the sII hydrate (cubic Fd3m space group) consists of 136 water molecules that form two types of cavities: the small, a pentagonal dodecahedron (512), and the large that is formed by © XXXX American Chemical Society

Received: January 15, 2013 Revised: May 2, 2013

A

dx.doi.org/10.1021/jp400449g | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

industrial communities.1 As a result, hydrates have been under consideration as a possible candidate material for storage and transportation of “energy-carrier” gases (e.g., methane,12 and hydrogen13,14) or the permanent storage of “greenhouse” gases (e.g., carbon dioxide15) inside sediments that are below the ocean floor. In addition, large amounts of methane hydrates have been discovered inside oceanic sediments and at permafrost regions, rendering them as possible future energy sources.16−18 The energy resources stored in the aforementioned hydrate deposits are estimated to be larger than the combined resources from all the other fossil fuels.19,20 Note, however, that the amount of methane stored in hydrates is an issue that is still under debate in the literature.21 During the recent years a significant effort has been dedicated to the study, both experimentally22−27 and theoretically/numerically,22,26−34 of methane production from permafrost and oceanic hydrate deposits. The aforementioned references are only some typical examples of the work that was reported in the particular research area. Providing a complete list of references is beyond the scope of the current study. Different production methods have been considered,17 including depressurization,29 thermal stimulation,28,30 and hydrate inhibitor injection.35 Essentially all production methods rely on the principle of shifting the prevailing reservoir conditions from stable conditions for hydrate to unstable conditions for hydrate, therefore leading to hydrate dissociation. Ultimately, the objective of all the proposed methods is to dissociate hydrates present within subsurface deposits in a controlled manner and produce the trapped methane gas in a safe, environmentally friendly, and economically viable way.36 Therefore, the gas saturation (i.e., the volume occupied by gas over the total pore volume) in a porous medium that results from the dissociation of methane hydrate is an important parameter. The ability to produce the gas phase is determined by whether the gas saturation is below or above a critical value for a particular system.37 For the case of primary oil production, the gas saturation corresponding to this critical limit is known as the “critical gas saturation” of the system.38−40 The resulting gas phase clusters (i.e., the gas bubbles that are distributed randomly in the porous medium reflecting the random initial distribution of the hydrate grains), in order to reach the production well, need to connect to form a “domain-spanning” gas cluster.41 Alternatively, the different gas clusters could get mobilized due to viscous/buoyancy forces.40 The amount of gas released during hydrate dissociation is equal to the amount of gas enclathrated (trapped) in the hydrate phase and depends on the temperature, pressure, and salinity conditions of the hydrate-bearing sediments. Therefore, it is essential to have methods to evaluate the gas saturation of different hydrate systems, based on the properties of the systems and the prevailing temperature and pressure conditions. The amount of hydrate in the system depends on the availability of the hydrate-forming methane gas inside the hydrate stability zone (HSZ),42 which is a domain with pressure such that the hydrate can be stable. Methane can be produced, either locally (e.g., biogenic origin) or far (e.g., thermogenic origin) from the HSZ, and subsequently be transported inside the HSZ by buoyancy forces.43,44 Depending on the source of methane, the porous domain can have different types of hydrate saturations. In particular, low hydrate saturations have been observed to result from biogenic methane, while high hydrate

saturations have been observed to result from thermogenic methane.42 Note, however, that this rule does not apply always. In particular, in Nankai Trough, offshore Japan, it has been found45 that methane hydrate fills the pore space, reaching saturations up to 80% despite the fact that methane gas is of microbial origin.46 The objective of the paper is to calculate the gas saturation resulting from hydrate dissociation in different types of domains/systems (i.e., nonporous and porous) with emphasis at the single pore and pore-network scales. Studies at the microscopic level can help to identify mechanisms occurring at the pore scale and guide in the effort to incorporate the observations to the continuum scale simulations, after performing the appropriate upscaling of the parameters involved. Experimental works examining hydrates at the pore-network scale were reported by a number of researchers.47−51 Hydrate formation inside a porous domain was studied numerically by Kang et al.52,53 using the lattice Boltzmann approach. Tsimpanogiannis and Lichtner54 extended the pore network approach that was also reported by Kang et al.53 and examined the effect of gravity gradients on hydrate formation inside pore networks. A pore-network-type simulator for hydrate dissociation was initially proposed by Tsimpanogiannis and Lichtner 55 and subsequently by others.56−60 Other hydrate studies at the pore-scale level involve the use of the discrete element method (DEM). Kreiter et al.61 used DEM in order to simulate hydrate growth in sediments with differing grain sizes. Jain and Juanes62 used DEM to simulate gas migration in sediments during the hydrate formation process. Bruggada et al.63 used DEM simulations to study the geomechanical behavior of hydrate saturated porous soils. Holtzman et al.64 developed a variational approach to simulate hydrates at the grain scale. This paper is organized as follows: Initially, we present the mathematical formulation for the case of a spherical solid hydrate particle in a bulk phase, namely, in the absence of interference with any solid material. Next, the case of a single hydrate particle inside a square, two-dimensional pore network is described. This is followed by the case of multiple hydrate particles that are distributed randomly in the pore network. We consider both cases of a network with all pore bodies/throats having equal size and a network with pore throats following an appropriate size distribution. Then, we discuss the comparison of the analytical results for all the aforementioned cases with simulations of hydrate dissociation inside pore networks, reported earlier.55 We also introduce two modifications in the pore network simulator that significantly improve the comparison with the analytical results. Finally, we end with the conclusions.

2. MATHEMATICAL FORMULATION: BULK PHASE Consider a spherical methane (CH4) hydrate particle with an initially known radius rH and volume VH, where the subscript H denotes the hydrate phase. The hydrate phase is placed inside an aqueous (H2O) liquid phase (denoted by the subscript l) of temperature T and pressure Pl. The hydrate/water system is located in a container of arbitrary shape that has a volume VD. Such a system is shown in Figure 1a for the particular case of the two-dimensional projection of a cubic container. We define the hydrate saturation, SH, of such a system as the ratio V SH ≡ H VD (1) B



Article

conditions of this study, methane forms a structure sI type hydrate.1 The sI methane hydrate consists of 46 hydrogenbonded water molecules that are grouped in two small (512) and six large (51262) cages. If all the eight cages of each hydrate unit cell are occupied by a single methane gas molecule (i.e., corresponding to stoichiometric conditions) results in a theoretical value for the hydration number nW = 5.75. However, experimental studies70−74 have reported hydration number for the sI methane hydrate that are nW > 5.75 (typically in the range 5.8−6.3). This result clearly indicates the existence of empty cages that are present in the stable methane hydrate. Similar conclusions were also reported from Monte Carlo molecular simulations.75 If we assume that the hydrate has density, ρH, and molecular weight, MWH, then, during the hydrate dissociation process, nH moles of hydrate (nH = VHρH/MWH) will produce ng moles of methane with nH = ng. If we also assume ng = λVg, we can calculate λ as follows: ρH ρH = λ= MWH MWM + nW MWW (6)

Figure 1. Schematic of a spherical hydrate particle in the bulk: (a) before dissociation and (b) after dissociation.

If the pressure Pl in the liquid phase is higher than the hydrate equilibrium pressure, Peq H , that corresponds to the particular temperature T, then the hydrate particle is thermodynamically stable. For hydrate stability to occur, we have also assumed that the liquid phase is saturated with dissolved methane, in order to prevent hydrate dissolution in the liquid phase that is driven by the unsaturated liquid phase.65−67 However, if Pl < Peq H , then the hydrate dissociates into liquid and gas phases with initial volumes Vl and Vg, respectively, such that VH = Vl + Vg (as shown in Figure 1b). The subscript g denotes the gas phase. Similarly, we can define the gas saturation of the system as

Sg ≡

where MWM = 16 is the molecular weight of methane and MWM = 18 is the molecular weight of water. By combining the previous eqs 3, 4, and 6, we can obtain the following cubic equation for the gas-bubble radius rg: A3rg 3 + A 2 rg 2 + A1rg + A 0 = 0

Vg VD

where the equation parameters Ai with i ∈ [0, 3] are given as follows:

(2)

The gas phase will create initially a spherical bubble with radius, rg, and volume, Vg [with Vg = (4/3)πrg3]. The pressure Pg in the gas bubble is given by the ideal gas law, corrected with the compressibility factor z, in order to account for the nonidealities of the gas phase at higher pressures. PgVg = zng RT

(3)

where R is the gas constant, T is the temperature, and ng is the number of gas moles. In the current study we used the component-specific equation of state (EoS) developed by Duan et al.68 to calculate the compressibility factor of methane. This particular EoS needs to be solved numerically and was selected because it gives very accurate results. The pressure inside the gas phase of the spherical bubble, Pg, is also related to the capillary pressure, Pc, and the pressure in the surrounding liquid phase, Pl, as follows: 2γ + Pl Pg = Pc + Pl = rg (4)

A3 =

4 πPl 3

(8a)

A2 =

8 πγ 3

(8b)

A1 = 0

(8c)

A 0 = −λzVHRT

(8d)

By solving analytically the cubic eq 7, we obtain the values for the gas-bubble radius rg, and therefore, we can compute the ratio of the gas-to-hydrate saturations as follows: Sg SH

=

⎛ rg ⎞3 =⎜ ⎟ VH ⎝ rH ⎠ Vg

(9)

This expression is one of the main results of interest to this study. This solution takes into account the effect of the capillary forces when the gas bubble radius is small.

where is γ is the interfacial tension between the gas and liquid phases, which is a function of temperature and pressure. In the current study γ is calculated using the methodology reported in ref 69. The dissociation of methane hydrate can be described through the following chemical reaction, which is valid for temperatures higher than 273.15 K, while for lower temperatures the water phase is in the form of ice after dissociation, instead of liquid phase: CH4 ·nW H 2O ⇔ CH4(g) + nW H 2O(l)

(7)

3. MATHEMATICAL FORMULATION: PORE NETWORK WITH IDENTICAL PORES In the previous section we examined the case of a spherical hydrate particle in the absence of any interactions between the hydrate grain and the solid walls of the container. In the current section we examine the case of the hydrate solid phase being inside a single pore body of a regular and uniform, pore network. Consider the schematic of such a simple pore network as shown in Figure 2. All pore bodies and pore throats are water-saturated, except the pore body that is located at the center of the pore network, which is hydrate-saturated. In order to obtain analytical solutions, we assume that all the M × N pores of the pore network have a spherical pore body

(5)

where nw is known as the hydration number (or hydrate number), which gives the number of water molecules per guest methane molecule. Under the temperature and pressure C



Article

Tsimpanogiannis and Lichtner55 extended the aforementioned definition of Scg to the case of hydrate dissociation in a pore network with pore bodies and throats obtained from appropriate size distributions. In the more general case of a pore network the different gas clusters (gas bubbles) get connected and form a “sample-spanning” gas cluster. The percolating gas cluster is significant because it facilitates the methane gas production, which is the primary production target during methane hydrate dissociation inside sediments. During gas production from hydrate deposits, the process should operate above the critical value, Scg. This production behavior is in contrast to the traditional case of primary oil production, where the process should operate below the critical value, Scg. This difference is a result of the fact that during primary oil production we are interested, mostly, in the liquid phase (and less in the gas phase). On the other hand, during gas production from hydrates we are interested only in the methane gas and not in the liquid phase (i.e., water).55,81 Consider a hydrate-saturated pore body that belongs to a pore network. If the radius of the pore-body, rp, is such that rp < rcp, then, upon hydrate dissociation, the gas bubble will remain isolated inside the pore body. However, if rp > rcp, then the gas bubble could grow outside the initial pore body and inside the neighboring pore bodies, provided that the minimum capillary threshold of the pore throats (connecting the pore bodies) will be surpassed. The pressure in the gas phase is now given by

Figure 2. Schematic of a N × N (N = 5) pore network with all pore bodies/throats having the same size.

with the same (known) radius rp (and volume Vp). In general M ≠ N; however, we used M = N in this study. Each pore body is connected with the neighboring pore bodies with horizontal and vertical pore throats that have cylindrical geometry and length Lt ≠ 0. The pore throats have radius rt. In this study we have also assumed that the radii of the pore body/throat are related through rp = βrt, with β > 1. Al-Raoush and Wilson76 reported values for β−1 in the range 0.314−0.787, corresponding to values of β in the range 3.18−1.27, respectively, for unconsolidated glass beads. For sandstones, values for β can range from, as low as, 1.92−2.75, that were reported from X-ray tomography,77,78 to as high as, 10, that were reported by Dullien and Dhawan.79 Even higher values, β ≈ 10−100, were reported for diatomites.80 In the case that β ≫ 1, the contribution of the pore throats to the total volumetric storage capacity (i.e., total volume, Vtotal PN ) of the pore network is very minimal and therefore can be ignored. Such is done for the remainder of this study. In the Supporting Information we present a brief discussion of possible corrections that are associated with contributions to the total volume of the pore network that originate from pore throats. Therefore, the main storage capacity of the pore network, in this study, is provided by the volume of the pore bodies: N total VPN =

Pg = Pct + Pl =

(11)

The aforementioned gas-growth process can be repeated many times (i.e., for as long as Pg ≥ Pct + Pl). The gas-phase growth will stop once the pressure in the gas phase becomes such that Pg ≤ Pct + Pl. At this point the gas−liquid interfaces of the gas clusters remain pinned inside the pore throats. Finally, at static equilibrium (i.e., when growth has stopped for all gas clusters) the gas phase would have a total volume given by

Vg = α*Vp

(12)

A value of α* < 1 indicates that the gas bubble is unable to expand (after hydrate dissociation) up to the point of reaching the pore-body walls. On the other hand, a value of α* > 1 indicates gas bubble expansion beyond the initial hydratesaturated pore body. By combining eqs 3, 6, 11, and 12, we obtain

M

∑ ∑ Vp(i , j) i=1 j=1

2γ 2γβ + Pl = + Pl rt rp

(10)

The pore throats provide the primary connections between neighboring pore bodies. Pore throats also provide the resistance to flow as a result of the capillary threshold that needs to be exceeded, Pct = (2γ/rt) cos ϑ (a zero contact angle, ϑ, is assumed in this study, which corresponds to water−wet solid surface), in order for the throats to be penetrated by an invading nonwetting fluid phase (e.g., methane gas). Tsimpanogiannis et al.81 defined the critical radius of a single pore body, rcp, that corresponds to the critical gas saturation, Scg, as the pore radius for which the gas phase, which results from the dissociation of hydrate, initially occupying the pore body, is enough to occupy completely the pore body. Namely, the produced-gas phase after hydrate dissociation expands and reaches the pore walls, expelling in the process all the liquid phase from the pore. The gas saturation at that point is identified as the critical gas saturation, Scg.

α* =

λzRTrp 2γβ + rpPl

Cf

(13)

where Cf is a correction factor that takes into account the contribution of pore throats. In the current study we consider that Cf = 1. In the Supporting Information we present a brief discussion for the case when Cf ≠ 1. Depending on the prevailing value for Pl, different values for α* can result. A low limiting case is obtained when Pl = Peq H and results in the minimum value of α*. Higher values of α* are obtained when Pl < Peq H . By recalling that VH = VP (i.e., the pore is initially completely hydrate-saturated), we can compute the ratio of the gas-to-hydrate saturations as follows: Sg SH D

=

Vg VH

= α* (14) dx.doi.org/10.1021/jp400449g | J. Phys. Chem. C XXXX, XXX, XXX−XXX


Article

In this part of the study we make the following two assumptions: (i) all pore bodies have equal pore radii (the same principle applies to all pore throats as well), and (ii) we ignore the contribution of pore throats to the total pore volume. By making the two aforementioned assumptions, we can obtain an analytical solution to the problem. Both assumptions can be relaxed. However, more complicated pore-network calculations are required in order to solve the problem. In real porous media, both pore bodies and pore throats can have sizes that follow appropriate size distributions, which can be measured experimentally.77,82

4. RESULTS AND DISCUSSION 4.1. Hydrate in Bulk Liquid Phase. We begin by examining the case of hydrate dissociation of a single, spherical, hydrate particle in the bulk aqueous phase (i.e., no effect of the solid walls of the container or of the solid matrix of a porous medium). Consider the schematic shown in Figure 3. The

Figure 4. Results for the dissociation of a spherical hydrate particle in the bulk. Ratio of gas-to-hydrate saturations, Sg/SH, as a function of temperature, T, for a series of pressures in the liquid phase.

higher gas expansions, and therefore, higher gas saturations in that spatial subdomain. For most cases (with the exception of higher temperatures and liquid pressures) the resulting gas saturation is higher than the initial hydrate saturation, namely Sg/SH > 1. 4.2. Hydrate in a Pore Network with Identical Pores/ Throats. 4.2.1. Case of a Single Hydrate-Occupied Pore. Next, we consider the case of a pore network with identically sized pore bodies, where only one of the pore bodies, chosen randomly, is occupied by methane hydrate. Figure 5 shows the ratio of gas-to-hydrate saturations, Sg/SH, as a function of the ratio of pressures, Pl/Peq H , and for

Figure 3. Schematic depicting the different zones of a hydratesaturated two-dimensional bulk domain.

schematic depicts a two-dimensional water-saturated bulk domain. Initially, a number of noninteracting, spherical, hydrate particles are placed in the system. Subsequently, a pressure gradient is applied to the system. The pressure in the liquid eq phase varies spatially and is given as Pl = εPeq H (TH ), where ε is a proportionality constant. The hydrate equilibrium pressure, Peq H 83 = f(Teq H ), is calculated using the correlation given by Moridis. Recall also that higher hydrate equilibrium temperature, Teq H, corresponds to higher hydrate equilibrium pressure, Peq H. Two pressure zones can be identified: (i) The hydrate stability zone (HSZ) is dominated by liquid pressures such that Pl/Peq H > 1, (i.e., ε > 1). Inside the HSZ the hydrate particle can remain stable, without any dissociation phenomena occurring, and (ii) the hydrate dissociation zone (HDZ), where hydrate dissociation occurs when a hydrate particle is placed in that zone. In this particular pressure zone 0 ≤ Pl/Peq H ≤ 1 (i.e., 0 ≤ ε ≤ 1). The two pressure zones are separated by the hydrate dissociation front, where ε = 1. We define as the hydrate eq dissociation front all the spatial domain such that Pl = Peq H (TH ). Figure 4 shows the ratio of gas-to-hydrate saturations, Sg/SH, as a function of temperature, T, for a series of pressures in the liquid phase, Pl. We consider only cases such that 0 ≤ ε ≤ 1.Values of ε > 1 are of no practical interest to this study, since the hydrate particles cannot dissociate. If we plot the ratio of gas-to-hydrate saturations, Sg/SH, as a function of the ratio of pressures, Pl/Peq H , and for various isotherms (Figure S1 in the Supporting Information), the following observations can be made. For a given isotherm, lower liquid pressures in a particular spatial subdomain result in

Figure 5. Results from a single hydrate grain inside a uniform pore network. Ratio of gas-to-hydrate saturations, Sg/SH, as a function of the ratio of pressures, Pl/Peq H , for 283.15 K and various pore body radii.

temperature equal to T = 283.15 K. The particular temperature was selected since it was used earlier in the pore network simulations reported in ref 55. We show results for a series of pore-body radii, rp, and β = 10. The same results as Figure 5 are also shown in the Supporting Information (Figure S2). However, the figure is focused on a more narrow pressure and saturation range and also using linear scale instead of logarithmic. We clearly observe that lower liquid pressures result in higher gas expansions and, therefore, higher gas E



Article

(that results from the hydrate dissociation of a single, hydrateoccupied, pore body) can occupy other pore bodies as well, in addition to the initial, hydrate-occupied one; and (ii) a subzone with liquid pressures such that Pl/Peq H → 1, which is located closer to the hydrate dissociation front, where the system is below the single-pore critical gas saturation. In this subzone the gas phase (that results from hydrate dissociation) is located at the center of the pore body as a spherical bubble and occupies the initial, hydrate-occupied, pore body partially. Figure 7 shows the ratio of gas-to-hydrate saturations, Sg/SH, as a function of temperature (for T > 273.15 K), along the eq methane hydrate equilibrium curve (Peq H , TH ). Results are shown for two different values of β (β = 5, 10). Additional cases of β (β = 2.5, 20) are reported in the Supporting Information (Figure S3). In all cases we observe that the solution obtained with eq 14 is identical to the bulk case, when the pore bodies have radii such that rp > 10−5 m. For the particular case of β = 2.5 this observation is valid for rp > 10−6 m. When we increase β, it results in shifting the Sg/SH curves to lower values. This effect becomes stronger for lower values of rp. For Teq H > 288 K and the larger values of rp, the saturation curves drop below the value of Sg/SH = 1, namely below the single-pore critical gas saturation. For the smaller values of rp, the temperature Teq H, where the gas saturation curve drops below the value of Sg/SH = 1, shifts to lower values. Again a stronger effect is observed for higher values of β. 4.2.2. Case of Multiple Hydrate-Occupied Pores. In the previous subsection we examined the case of a single hydrateoccupied pore body inside a pore network that undergoes hydrate dissociation. The pore body can be anywhere between the hydrate dissociation front and the point where the pressure, Pl, has the minimum value. Depending on the prevailing value for Pl, different values for α* could result. Two limiting cases can be obtained: (i) a minimum value of α* is obtained when Pl = Peq H (i.e., on the hydrate dissociation front), and (ii) a maximum value of α* is obtained when Pl = 0. In a more general case, multiple hydrate clusters are present inside the pore network. In addition different values of pressure, Pl, exist at the different positions within the pore network. In that case, eq 14 can be used to calculate analytically an average value for the ratio of gas-to-hydrate saturations, ⟨Sg/SH⟩analytical, for the entire pore network as follows:

saturations. Notice that for the particular temperature examined and for cases of pore bodies with approximate radii such that rp > 10−4 m and pore throats with radii rt > 10−5 m the solution obtained with eq 14 is identical to the case of bulk phase [given by eq 9]. We observe from Figure 5 that for pore-body radii such that rp > rcp (for the particular hydrate equilibrium conditions examined: rcp = 1.688 × 10−7 m) the ratio of gas-to-hydrate saturations is Sg/SH > 1, for all values of Pl/Peq H . This clearly indicates that the hydrate-occupied pore is going to be above the single-pore critical gas saturation, Sgc, when hydrate dissociation occurs. The observation is valid, no matter where the hydrate-occupied pore is located within the pore network, as long as it is within the hydrate dissociation zone (see also the schematic of Figure 6a). On the other hand, for pore-body radii such that rp < rcp two clearly distinct subzones can be identified within the hydrate dissociation zone, as can be shown in the schematic of Figure 6b: (i) a subzone with lower liquid pressures which is located closer to the low-pressure end, where the system is above the single-pore critical gas saturation. In this subzone the gas phase

analytical

⟨Sg /SH⟩

=

n ∑i = 1 αi*

n

(15)

If the pressure within the pore network varies in the range [Plow, Peq H ], then n denotes the number of points that are considered for the summation. If n is taken large enough, the analytical calculation converges to a constant value. Figure 8 shows the ratio of gas-to-hydrate average saturations values calculated (using eq 15) as a function of ΔPD = 1 − Plow/ Peq H , and for three different pore-body radii, at 283.15 K. Shown also are the maximum and minimum values discussed previously. The calculation is performed using the following steps: • Define a value for Plow, for the low-pressure-end of the pore network. Use Peq H for the high-pressure end of the pore network. Therefore, define a value for ΔPD, applied at the two opposite ends of the pore network. • Calculate the maximum value for α* using eq 13 by taking Pl = Plow.

Figure 6. Schematic of a porous medium with a hydrate stability zone (HSZ) and a zone where the hydrate has dissociated. (a) The entire hydrate dissociation zone (HDZ) is above the single-pore critical gas saturation. (b) The HDZ is divided in two subzones, above and below the single-pore critical gas saturation. F



Article

Figure 7. Results from a single hydrate grain inside a uniform pore network. Ratio of gas-to-hydrate saturations, Sg/SH, as a function of temperature (along the hydrate equilibrium curve) and various pore radii. Shown also with circles are the results from the bulk case. The horizontal dashed lines denote the single-pore critical gas saturation.

Figure 9. Ratio of gas-to-hydrate saturations, ⟨Sg/SH⟩, as a function of the dimensionless pressure drop across a 100 × 100 pore network, ΔPD, at 283.15 K, and for three different pore-body radii. Comparison of the analytical solutions (solid lines) and the results from the refined pore network simulator (triangles).

Figure 8. Analytical solution results: Ratio of gas-to-hydrate overall saturations, ⟨Sg/SH⟩, as a function of the dimensionless pressure drop across a pore network, ΔPD, at 283.15 K and for three different porebody radii. For each pore body radii the lower dashed line corresponds to a minimum value of ⟨Sg/SH⟩ obtained when Pl = Peq H , the upper dashed line corresponds a maximum value of ⟨Sg/SH⟩ obtained when Pl = 0, and the solid line corresponds to the average value for the pore network.

realizations. Very good agreement is observed between the two different approaches that were considered. The agreement becomes even more noteworthy, if we take into account the simplicity of the analytical calculation when compared with the complexity of the pore network simulator. This issue is further discussed in the following section. 4.3. Hydrate in a Pore Network with Pore Bodies/ Throats Following Size Distributions. Tsimpanogiannis and Lichtner55 developed a pore network simulator for hydrate dissociation and reported results mostly for the case of low hydrate saturations in pore networks with low permeability. Such conditions correspond to the case of a class 4 type of natural hydrate deposit (i.e., with hydrate saturation less than 5%, but locally up to 20%), which is the common type of oceanic hydrate deposit.84,85 The simulator assumes instantaneous dissociation (i.e., thermodynamic equilibrium model) of the hydrate phase, once the pressure drops below the hydrate equilibrium value, Peq H (T), for a given temperature T. This could be a realistic assumption in cases when small hydrate clusters, namely each one having a size of a few pore bodies, are distributed throughout the entire porous domain. However, this

• Calculate the minimum value for α* using eq 13 by taking Pl = Peq H. • Plot α* vs Pl, using eq 13. The pressure, Pl, within the pore network varies in the range [Plow, Peq H ]. • Calculate ⟨Sg/SH⟩analytical for the particular pressure range using eq 15. The calculated value corresponds to the analytical solution for the ratio of gas-to-hydrate saturations that would result for the given pore network when a ΔPD across the network is applied having a minimum pressure equal to Plow. Notice that the minimum value for α* is constant for all values of Plow considered, while the maximum value for α* is variable and depends on the value of Plow selected. Figure 9 shows the comparison between the ratio of gas-tohydrate average saturations calculated using the analytical solution from eq 15 and the results obtained from the detailed pore network simulator. For the pore network results we used two-dimensional networks of size 100 × 100. The results for every point shown on the figure are averaged over 50 G



Article

were subscript k denotes the phase under consideration [k ∈ (g, H), for gas and hydrate, respectively], and ωk is a phase function that is defined as follows:

assumption needs further investigation for cases the hydrate clusters become significantly larger, and therefore, instantaneous dissociation may be unrealistic to describe the physics of the problem. Examining the case of higher hydrate saturations could require possibly the inclusion of hydrate dissociation kinetics86 into the simulator, where the driving force for the hydrate dissociation is proportional to the difference between the equilibrium fugacity of hydrate at a given temperature and the gas fugacity. Such an effort is currently in progress in order to evaluate the effect of dissociation kinetics on the resulting gas-phase growth patterns and the various saturations at the pore network scale. Reported results from field-scale numerical modeling studies87 indicate that the equilibrium and kinetic models for methane dissociation in a porous medium are almost indistinguishable. On the other hand, the hydrate dissociation kinetics appear to be important for short-term processes and core-scale simulations.87−89 Additional discussion on the use of nonequilibrium models is presented by Vafaei et al.90 In the current study the basic aspects of the pore-network simulator remain unchanged, and thus no further discussion is offered here, since the details are given in ref 55. In this study, however, the amount of gas that results from the dissociation process is calculated with significantly higher accuracy than in ref 55. Here, we consider the detailed thermodynamics and solve the hydrate dissociation problem for a given temperature and pressure in order to calculate the hydrate cage occupancies and, therefore, the exact amount of gas that is released during the dissociation process. Please recall that clathrate hydrates are nonstoichiometric compounds. Namely, all hydrate cages are not required to be fully occupied by the guest molecules in order for the hydrate structure to be stable. It is well-known from reported experiments70−74 that for the case of methane hydrates the cage occupancies, θi, corresponding to the small and large cages are θi < 1 (with i: S and L denoting the small and large hydrate cages, respectively). In this study the initial hydrate saturation is distributed randomly in the pore network after the value for the initial hydrate saturation is decided. A similar approach was followed in ref 55. Therefore, it is essential to repeat the simulations using different realizations for each one of the hydrate saturations examined. We consider isothermal simulation conditions that are similar to those reported in our earlier study.55 In particular, we examine methane hydrate dissociation at 283.15 and 275.15 K. We use square, two-dimensional, pore networks with similar sizes as well as similar pore body/throat size distributions with β = 10. The following two groups of parameter values are used 5 eq during the simulations: (i) Teq H = 283.15 K, PH = 72.43 × 10 −3 3 Pa, z = 0.8577, γ = 65.79 × 10 N/m, λ = 7.4487 × 10 ; (ii) eq 5 Teq H = 275.15 K, PH = 32.00 × 10 Pa, z = 0.9282, γ = 68.80 × −3 3 10 N/m, λ = 7.3395 × 10 . 4.3.1. Average Overall (Total) Saturations. Once the PN simulation is completed for each separate realization, overall gas saturations are calculated using the entire pore network. Both gas and hydrate saturations are defined as the volume fractions of gas and hydrate, respectively, by considering the volume of the entire domain of the pore network as follows: Sk =

N N ∑ j = 1 ∑i = 1 ωk(i , j)Vk(i , N N ∑ j = 1 ∑i = 1 Vp(i , j)

⎧1 if pore‐body (i , j) contains phase k ω k(i , j) = ⎨ ⎩ 0 otherwise

(17)

In a final step, the overall saturations for the phase k are averaged over the different realizations in order to obtain the average overall saturation for the particular phase, ⟨Sk⟩. Figure 10 shows the ratio of gas-to-hydrate average overall saturations, ⟨Sg/SH⟩, as a function of the dimensionless pressure

Figure 10. Ratio of gas-to-hydrate average overall saturations, ⟨Sg/SH⟩, as a function of the dimensionless pressure drop across the pore network, ΔPD, at 283.15 K (red) and 275.15 K (blue) and for various initial overall average hydrate saturations (⟨SH⟩ = 1%, 3%, 6%, 9%, and 12%). Solid lines denote the results from the analytical solution. Results are obtained from 100 × 100 pore networks and are averaged over 50 realizations of randomly distributed hydrate of given values of ⟨SH⟩ (rp = 5 × 10−7 m).

drop across the pore network, ΔPD, at temperatures 283.15 and 275.15 K and for various initial average overall hydrate saturations (⟨SH⟩ = 1%, 3%, 6%, 9%, and 12%). The dimensionless pressure drop across the pore network, ΔPD, is defined as follows (see also Figure 6 for a schematic): Phigh − Plow ΔPD = Phigh (18) In principle, Phigh could have values higher than Peq H . Such approach was followed in ref 55. In the current study, however, we report results that were obtained by using Phigh = Peq H . In this case, all the hydrate phase in the pore network dissociates during the simulations. The results presented in Figure 10 are obtained from simulations in 100 × 100 pore networks and are averaged over 50 realizations. The average pore body size is equal to rp̅ = 5 × 10−7 m and the pore throats have radii in the range rt ∈ [5.01 × 10−8 − 4.99 × 10−8] m. As it should be expected, the results obtained from the different initial hydrate saturations collapse to a single curve when plotted as ⟨Sg /SH⟩ vs ΔPD

j)

(19)

Deviations from the collapse on the single curve are very small. In particular, for the case of 283.15 K the maximum deviation is less than 3% for all ΔPD < 0.99, it is less than 1% for all ΔPD

Gas Saturation Resulting from Methane Hydrate Dissociation in a

Recommend Documents