Sensitivity Analysis of Parameters Governing the Recovery of Methane from Natural Gas Hydrate Reservoirs

Naturally occurring gas hydrates are regarded as an important future source of energy and considerable efforts are currently being invested to develop methods for an economically viable recovery of this resource. The recovery of natural gas from gas hydrate deposits has been studied by a number of researchers. Depressurization of the reservoir is seen as a favorable method because of its relatively low energy requirements. While lowering the pressure in the production well seems to be a straight forward approach to destabilize methane hydrates, the intrinsic kinetics of CH 4-hydrate decomposition and fluid flow lead to complex processes of mass and heat transfer within the deposit. In order to develop a better understanding of the processes and conditions governing the production of methane from methane hydrates it is necessary to study the sensitivity of gas production to the effects of factors such as pressure, temperature, thermal conductivity, permeability, porosity on methane recovery from naturally occurring gas hydrates. A simplified model is the base for an ensemble of reservoir simulations to study which parameters govern productivity and how these factors might interact.


Introduction
Naturally occurring gas hydrates are regarded as an important future source of energy and considerable efforts are currently being invested to develop methods for an economically viable recovery of this resource.The recovery of natural gas from gas hydrate deposits has been studied by a number of researchers [1][2][3][4][5][6].Depressurization is seen as a favourable method because of its lower energy requirements [7].While lowering the pressure in the production well seems to be a straight forward approach to destabilize methane hydrates, the intrinsic kinetics of CH 4 -hydrate decomposition and fluid flow lead to complex processes of mass and heat transfer within the deposit.In order to develop a better understanding of the processes and conditions governing the production of methane from methane hydrates it is necessary to study the sensitivity of gas production to the effects of factors such as pressure, temperature, thermal conductivity, permeability, porosity on gas production.The main purpose of this paper is therefore to carry out a sensitivity study on methane recovery from naturally occurring gas hydrates using a reservoir simulation in order to determine which parameters govern productivity and how these factors might interact.In this model, the three primary mechanisms involved in the gas hydrate decomposition are considered to be: (1) heat transfer; (2) gas-water fluid flow; and (3) kinetics of hydrate decomposition.

Sensitivity Analysis-Theoretical Background
With the development of new computational techniques and robust numerical calculations in the last decades, reservoir simulation has evolved into a cost-effective tool for industries, not only for testing new production designs but also for maximizing the production by optimizing the process variables.It is in this context that sensitivity analysis can contribute to phenomenological understanding, key factor determination, and assessing interactions between factors, process optimization and production forecasting [8].Sensitivity analysis is a technique for systematically changing factors in a model to determine the effect of such changes in either one or several response variables [9] and has been used in several numerical reservoir studies [8,10,11].Two approaches were selected to carry out the sensitivity analysis: (1) one factor at a time (OFAT); and (2) factorial design.In this paper, the purpose of the sensitivity analysis will be the determination of the relative importance of selected process variables and physical properties on the cumulative gas produced from a hypothetical gas hydrate deposit.

OFAT
OFAT his is one of the most common methods to carry out sensitivity analysis by selecting a baseline for each factor and then sequentially varying each factor over its range with the other factors fixed at the initial (baseline) level [12].Based on this approach any change observed in the output will unambiguously be due to the single factor changed.This increases the comparability of the results and minimizes the risk of a numerical model not converging, which is more likely to happen when several input factors are changed simultaneously [9].The major disadvantage of this strategy is that it does not consider any interaction between the factors.An interaction, in this case, is characterized by the failure of one factor to produce the same effect on the response variable when another factor is also changed [12].

Factorial Design
In factorial design, all factors are varied together taking all possible combinations in the different levels of each factor.The analysis of the results is usually carried out using analysis of variance (ANOVA); which is a collection of statistical models in which the observed variance in a particular variable is partitioned into components attributable to different sources of variation [12].The main advantage of this approach, when compared to OFAT, is that interaction effects between variables are considered.It allows the effect of a factor to be estimated at several levels of other factors and is thus more powerful than OFAT.The purpose of ANOVA is to test differences in means (for variables or groups) for statistical significance.This is achieved by partitioning the total variance in a measured outcome into its sources: The components that are due to differences between means (SS Effects); which means variance that can be explained, such as by a regression model or an experimental treatment assignment and a component that is due to true random error (SS Error), which means a variance that cannot be explained.These latter variance components are then tested for statistical significance using the F-test or the p-value [12].In this paper, all the outcomes come from a simulation, thus, no true random error exists which causes ANOVA analysis to fail.In order to overcome this problem, it is assumed that high order interactions are negligible and combine their mean squares to estimate the error.This approach is based on the "sparsity of the effects principle" [12] which states that most systems are dominated by some of the main effects and lower interactions and higher order interactions are negligible.The calculation procedure will be presented in the subsequent sections.

Modeling CH 4 Recovery from CH 4 -Hydrate Reservoirs
To study the characteristics of CH 4 -hydrate decomposition in porous media we developed a model based on the following considerations: Gas hydrate decomposition can be represented by the following kinetic reaction: CMG STARS™ [13] (Computer Modelling Group-Advanced Processes & Thermal Reservoir Simulator, Calgary, AB, Canada) was chosen as the tool to solve the mass and energy balances.Its suitability to represent methane production from gas hydrates in porous media and has been examined by other authors who validated it against other comparable software [ [14][15][16][17].

Parameter
Chen et al.

Reservoir Parameters and Component Properties
In the simulation of gas hydrate several parameters are considered as being important for the production of methane from CH 4 -hydrate deposits.The most important are: geometry of the reservoir (dimensions and shape), reservoir properties (porosity, permeability, etc.), initial conditions (temperature, pressure, gas/water saturation, and hydrate saturation), operational conditions (bottom hole pressure (BHP) and well radius).Table 1 summarizes the available data for methane production from gas hydrate reservoirs.Values from the literature were selected to create a base case and subsequently use value limits to verify their impact on methane recovery using a sensitivity analysis.The data selected for the base case and their variations will be presented in the subsequent sections.
As in the case for reservoir properties, a literature review was conducted for gas hydrate physical properties and it is summarized in Table 2 [1,5,16,[18][19][20]. From those values, we selected representative data to be used in the base case.

Governing Equations and Model Design
The model was designed as a radially symmetrical model of a gas hydrate deposit.It includes all three rate-controlling mechanisms (multiphase fluid flow in porous media, kinetics of decomposition and heat transfer) that govern CH 4 production from naturally occurring gas hydrates [16]: Mass balance for CH 4 : ( Mass balance for water: Mass balance for CH 4 -Hydrate: (3) Fluid flow velocities are described by the following equations: The velocity of gas along r-direction: The velocity of gas along z-direction: The velocity of water along r-direction: The velocity of water along z-direction: (7) The rates of gas and water production and the rate of hydrate decomposition are represented by and , and are calculated using the Kim-Bishnoi model [25]: The energy balance is described by the following equations: (11) where: Energies 2014, 7 2154

Gas Hydrate Decomposition in STARS™
As mentioned before, the kinetics of methane hydrate decomposition follows the Kim-Bishnoi model [25].The value for the heat of decomposition is taken from Liu et al. [1]: ; ΔH R = 54.7 kJ/mol (13) However, this expression needs to be rearranged to conform to the notation used in STARS™ [15,16]. ; Kinetic parameters (, , E) were determined experimentally by Clarke and Bishnoi [26].For pure components the fugacity of methane can be calculated as: f CH4 = φP; f e = φP e .Since STARS™ does not allow fugacity calculation the system is assumed to behave as an ideal gas, and we thus assigned the value φ = 1: The equilibrium pressure equation is determined by regression from experimental data compiled by Sloan and Koh [27] (see Figure 1): (17) and: (18) (19) Thus: (20) Gas hydrate decomposition affects reservoir properties such as permeability and porosity since solid gas hydrate partially blocks the reservoir pore space, thus decreasing its permeability.During the decomposition process the pore space is freed from solid gas hydrate, thus increasing its permeability, which favors fluid transport though the reservoir.To take into account the effect of hydrate saturation in the pore space, STARS™ offers the option to determine the reservoir absolute permeability, k 0 , with each time step using a Carmen-Kozeny type formula: where ϕ is the effective fluid porosity; and ε is an empirical parameter.On the other hand, the relative permeability and capillary pressure were obtained from Hong [17] who modified the equations of Van Genuchten [28] and Parker et al. [29] to incorporate the presence of a hydrate phase.[27] to determine the equilibrium pressure P* (Equation ( 19)).
Gas and aqueous relative permeability are thus described by the following equations: (22) (23) Capillary pressure between gaseous and aqueous phases is expressed as: (24) where:

Geometrical Reservoir Characteristics and Exploitation Period
Following the simulation carried out by Hong [17], the chosen case study was a reservoir of cylindrical shape with a single production well located at the center.The pore space is occupied by mobile phase (gas and water) and immobile phase (methane hydrate).The reservoir formation is considered to be a homogeneous and isotropic medium.In the numerical simulation the system is allowed to lose or gain heat from the surroundings in order to make the model realistic.A schematic representation of the system is presented in Figure 2. A trial and error method was used to determine the time of the simulation in order to ensure the development of the decomposition profile that makes best use of the available hydrates.Three periods of time were evaluated: 8, 12 and 20 years.The results of these simulations are shown in Figures 3-5.

Results and Discussion
Before starting analyzing the different results, it is important to highlight that the purpose of this paper is to carry out a systematic sensitivity analysis of variables, thus, the selection of parameter's values is not done based on a single type of formation encountered on Earth, but rather, a collection of multiple combinations that although might create "non-possible" hydrate reservoirs will also contain all the hydrate reservoirs found in the different regions, allowing the authors to unveil the what the most important variables are and its different interactions among each other.

Base Case
The literature review presented in previous sections provided the parameter values that were selected to create the base case for a sensitivity study to investigate the influence of the chosen parameters on gas recovery from methane hydrates (Table 4).
Figure 6 presents the cumulative amount of recovered methane and the rate of gas production obtained from the simulation of the base case.In this scenario, the cumulative amount of gas recovered from the reservoir increases almost constantly during the 12 years of operation.However, the rate of gas production reaches its peak in the first years of operation and subsequently declines.Figures 7 and 8 show the temperature and CH 4 -hydrate concentration distribution in the reservoir base case at the end of the operation.As noted, gas hydrates have decomposed in the area close to the over-and under-burden.Due to the endothermic nature of the hydrate decomposition process the temperatures in the reservoir decrease, slowing the decomposition process, while at the same time heat is provided from the surrounding over-and under-burden, which in turn accelerates the decomposition process.To ensure that the results of the models are comparable through all modifications of the reservoir parameters we kept the total molar amount of gas hydrate in the reservoir constant.Two parameters-porosity and thickness of the hydrate layer-deserve special attention since any modification in these values causes a change in the formation pore volume and volume of the reservoir.In these two cases, it was necessary to modify the geometry of the system to keep the total volume of the reservoir constant.

OFAT
Table 5 summarizes the results of the OFAT.The parameters chosen for this sensitivity study were systematically changed to extreme values and the results compared with the base case.
Table 5. Results of the one factor at a time (OFAT) analysis.
Radius is increased by a factor of 2 Figures 9-15 present the cumulative recovered methane and gas rate production obtained from the simulations of the variables when changed one at a time.As can be noted from these results, the main variables that produce significant changes in the total volume of gas that can be produced in a 12-year production period are absolute permeability (κ) and BHP.
Figures 10 and 12 show how changes in these parameters result in distinctive changes in gas rate production.In the case of the absolute permeability, the gas rate production curve increases very sharply in the first year reaching a maximum and subsequently a quick drop in the gas rate; then during the following years, the gas rate seems to reach a plateau.BHP had a negative effect on the cumulative gas production of methane, since there is less driving force to produce gas.The gas rate production is delayed, producing a maximum peak at approximately five years to the subsequently decrease.
Variables such as thermal conductivity of the rock, initial pressure and well radius had a smaller effect on the response variable (cumulative gas production) indicating that their significances are low in comparison to absolute permeability and BHP.Gas rate production curves (Figures 9,11 and 15) showed similar trends to the ones presented by the base case.
As noted above, changes in porosity and thickness of the hydrate layer affect the formation pore volume and the volume of the reservoir, respectively.In the case of a decreased porosity it was necessary to increase the thickness of the hydrate layer by a factor of 2.5 to compensate the change in porosity from 0.5 to 0.2.In the case of a reduced thickness of the gas hydrate it was necessary to increase the radius of the reservoir by a factor of 2.5 to compensate the change in thickness from 14 m to 7 m.    Figure 13 shows that in the case of porosity the cumulative recovered methane is higher than in the base case.Although, common sense would suggest that low porosity will produce less gas, the results show that in this case the cumulative production is higher.This case shows an increased ratio of formation rock against hydrates that results in a higher heat retention of the formation rock that subsequently result in a higher temperature in the host formation and subsequently temperatures remain higher due to the heat capacity of the rock, which in turn stimulates even more hydrate dissociation.In the case of reducing the thickness of the hydrate layer, the radius of the reservoir was increased in order to keep the same initial volume.The result of this change was less cumulative methane production and is shown in Figure 14.In this case gas from hydrate decomposition had to travel a longer distance through the host rock and an increased time was needed for the pressure change to propagate through the formation.

Full Factorial Design
For our full factorial design we chose six parameters: BHP, initial pressure, thermal conductivity of the rock, absolute permeability, irreducible water saturation and well radius.The effects of each factor with respect to the response variable were studied in an ANOVA.The design bases for this analysis are: • Six factors and two levels; this design constitutes a total of 64 simulations (2 6 ); • Response variable: cumulative production of CH 4 ; • Level of significance (α): 0.05; • Replicates: one (number of results per run).
Table 6 summarizes the factors and levels that are used in this analysis.Table 7 summarizes the results of the 64 simulations required to study this system.A value of −1 represents a low level of the factor whereas a value of +1 represents a high level.

Factors
Run label
The chosen parameter values were BHP = 4000 kPa; initial pressure = 6000 kPa; thermal conductivity = 0.5 W/(m•K); absolute permeability = 300 mD; irreducible water saturation = 0.2; well radius = 0.086 m.This run is labeled "ad" because the factors a (BHP) and d (absolute permeability) are at the high level.A similar structure is followed by Runs 20 and 37; where the highest levels are respectively: absolute permeability and well radius; thermal conductivity, absolute permeability and irreducible water saturation.
As mentioned before, this system is an unreplicated 2 6 factorial design.With only one replicate there is no internal estimate of error, which causes ANOVA analysis to fail.To overcome this problem, it is assumed that high order interactions are negligible (sparsity of effects principle) and combine their mean squares to estimate the error [12].
To carry out this type of ANOVA it is necessary to calculate its main components: Sum of squares (Sum Sq.) due to each source represent the difference between the various levels of each source and the grand mean; degrees of freedom (d.f.) represents the degree of freedom associated to the sum of squares; mean squares (Mean Sq.) of each source represent the variance of each source and is obtained by dividing the sum of squares by the associated d.f.F-test is obtained from the ratio between variance of the source and the variance of the error.P-value is the statistical test to determine if a factor is statistically important, in this case, P-value is less than the level of significance (α).Table 9 presents the results obtained from the ANOVA analysis.
From Table 9, using the 0.05 level of significance, it is possible to deduce that the factors BHP, K (thermal conductivity of the rock) and κ (absolute permeability), and the interactions BHP-K (BHP-thermal conductivity), BHP-κ (BHP-absolute permeability) and K-κ (thermal conductivity-absolute permeability), all have a value of p < 0.05 and thus make statistically important contributions to the response variable.
It is important now to rule out higher order interactions.Montgomery [12] proposes a method of analysis which is based on the normal probability plot of the estimates of the effects.In this method, the effects that are negligible are normally distributed, with mean zero and will tend to fall along a straight line on this plot, whereas significant effects will have non-zero means and will not lie along the straight line.In order to calculate the effects it is necessary to determine the contrasts for each effect.The general method is as follows In general to estimate the contrast for effect the AB…K, expand the right hand side of the following equation: (28) important (parallel lines) since the differences between the well radii are consistent for each value of BHP; however, the interaction BHP-κ is statistically important (non-parallel lines) since at 0.1 mD the average cumulative gas produced is almost the same at the two BHP levels, but at 300 mD there is a large difference.An additional interpretation of the effects is possible; since the effects of initial pressure, irreducible water saturation and well radius and all their interactions are negligible it is possible to discard these variables so that the design becomes 23 with 8 replicates, this method is known in the literature as hidden replication.The ANOVA using this simplification is presented in Table 10.The results of this analysis confirm the conclusions obtained above.

Conclusions
The OFAT sensitivity study and the full factorial design analyzed using ANOVA showed that reservoir absolute permeability (κ), BHP and the thermal conductivity of the rock (K) had the most significant effects on the recovery of methane from gas hydrates.In addition, the full factorial design w allowed to established important interactions such as BHP-K, κ-K, BHP-κ and κ-K-BHP that contribute to the response variable which cannot be concluded from the OFTA.
Considering the levels, by which these factors contribute to productivity, the most significant factors, by far, are those controlling the transport of gas to the well (κ and BHP).In comparison to factors governing the rate of transport, factors governing the rate of hydrate decomposition are of less importance, except for heat transfer (K) from the country rock to the hydrate deposit.
The design of this study will help to estimate the upper limits of productivity in natural gas hydrate deposits.In the field, further limitations will apply, such as trapped gas bubbles, inhomogeneous hydrate distribution or anisotropic permeability.

Figure 2 .
Figure 2. Schematic representation of the hydrate reservoir in Computer Modelling Group-Advanced Processes & Thermal Reservoir Simulator (CMG STARS).

Figure 3 .
Figure 3. Concentration of solid methane hydrate in a model reservoir after eight years of production.Most methane hydrates still remain in the reservoir.

Figure 4 .Figure 5 .
Figure 4. Concentration of solid methane hydrate in a model reservoir after twelve years of production.Significant amounts of methane hydrate remain in the reservoir.

Figure 6 .
Figure 6.Cumulative gas production in the base case scenario.

Figure 7 .
Figure 7. Temperature distribution in the reservoir base case at the end of the operation.

Figure 8 .
Figure 8. CH 4 -hydrate concentration distribution in the reservoir base case at the end of the operation.

Figure 9 .
Figure 9. Influence of thermal conductivity of the rock on gas production.

Figure 10 .
Figure 10.The influence of absolute permeability on gas production.

Figure 11 .
Figure 11.Influence of initial pressure on gas production.

Figure 12 .
Figure 12.Influence of bottom-hole pressure (BHP) on gas production.

Figure 13 .
Figure 13.Influence of porosity on gas production.

Figure 14 .
Figure 14.Influence of thickness of the hydrate layer on gas production.

Figure 15 .
Figure 15.Influence of well radius on gas production.

•
Reservoir is an unsteady-state open system; • Phases involved are: gas, aqueous, hydrate; • Components: CH 4 , H 2 O, CH 4 -hydrate; • Darcy's law will be used to describe the fluid flow in porous media; • Reservoir is a homogeneous and isotropic medium; • Gas hydrate molecule is assumed to be CH 4 •5.75H 2 O.

Table 2 .
Literature review of gas hydrate properties.

Table 3 .
Physical properties of the compounds involved in the formation of methane hydrates.

Table 4 .
Data for gas hydrate reservoir base case.

Table 6 .
Factors and levels of the factorial design.

Table 8 .
Label structure for simulation runs.

Table 10 .
ANOVA results-23 factorial design "hidden replication". : H Heat of hydrate decomposition unit bulk volume (J/m 3 •s) Q in Direct heat input bulk volume (J/m 3 •s)