Numerical Simulation of the Depressurization Process of a Natural Gas Hydrate Reservoir : An Attempt at Optimization of Field Operational Factors with Multiple Wells in a Real 3 D Geological Model

Zhixue Sun 1, Ying Xin 1, Qiang Sun 1,*, Ruolong Ma 2,*, Jianguang Zhang 1, Shuhuan Lv 1, Mingyu Cai 1 and Haoxuan Wang 1 1 School of Petroleum Engineering, China University of Petroleum (East China), Qingdao 266580, China; upcszx@upc.edu.cn (Z.S.); m15764220169@163.com (Y.X.); Eduzjg@163.com (J.Z.); 15553243290@163.com (S.L.); cmy9724@163.com (M.C.); 13506428396@163.com (H.W.) 2 College of Energy Resources, Chengdu University of Technology, Chengdu 610059, China * Correspondence: sun-q001@163.com (Q.S.); mrl0405@126.com (R.M.); Tel.: +86-532-8698-1170 (Q.S.); +86-28-8407-9071 (R.M.)


Background
The global temperature of the Earth's atmosphere is rising due to human emissions of greenhouse gases, especially carbon dioxide produced by fossil fuel combustion [1].Since the beginning of the Western Industrial Revolution, the acquisition of clean, secure and sustainable energy sources has been the future and foundation of the world's increasing economic growth [2].As a substitute for Energies 2016, 9, 714 2 of 20 fossil fuel and potential source of energy, natural gas hydrates are drawing global scientists' attention and playing an important part in resolving the energy crisis and climate change issues [3]. Figure 1 indicates that they are wildly distributed around the global, being found in particularly large quantities as sediments deposited on the deep ocean floor and as permafrost on the continental plateaus [4].Natural gas hydrates are crystalline solids composed of small gas molecules and water in which the gas molecules are extremely compressed, making them denser than typical fluid hydrocarbons [5].The rough estimated volume of gas that could be released from the hydrate reservoirs all over the world already exceeds that of known traditional gas resources [6].It is probable that we will use the new environmentally friendly fuel source to meet the increasing demand for energy in the next two decades.In terms of various areas worldwide, Figure 1 compares the reserve volume of methane hydrate.Consequently, we can draw the conclusion gas hydrates are one kind of new energy with great promise in the future [1].
Energies 2016, 9, 714 2 of 20 and playing an important part in resolving the energy crisis and climate change issues [3]. Figure 1 indicates that they are wildly distributed around the global, being found in particularly large quantities as sediments deposited on the deep ocean floor and as permafrost on the continental plateaus [4].Natural gas hydrates are crystalline solids composed of small gas molecules and water in which the gas molecules are extremely compressed, making them denser than typical fluid hydrocarbons [5].The rough estimated volume of gas that could be released from the hydrate reservoirs all over the world already exceeds that of known traditional gas resources [6].It is probable that we will use the new environmentally friendly fuel source to meet the increasing demand for energy in the next two decades.In terms of various areas worldwide, Figure 1 compares the reserve volume of methane hydrate.Consequently, we can draw the conclusion gas hydrates are one kind of new energy with great promise in the future [1].
Figure 1.Gas hydrate resource potential by global region, adapted from Johnson [7], with permission from GRID-Arendal, 2015.
The following equation describes the essence of hydrate formation and dissociation and the direction of the reaction: where NH is the hydration number and NH = 5.75 is the complete hydration number, while NH = 6.0 corresponds to the average number.ΔH is the enthalpy of the formation and dissociation reactions [6].
The hydration process in Equation (1) has also been thoroughly studied in the past few years as a promising method to separate methane from gaseous mixtures [8][9][10][11].Recently, several ways of developing hydrate reservoirs (Figure 2) have been put forward, including thermal stimulation [12][13][14], injection of hot brine, steam or hot water to heat the hydrate reservoir to exceed the dissociation temperature; depressurization [15][16][17], to lower the original pressure of hydrate reservoir under the dissociation pressure at a specific temperature; thermodynamic inhibitor injection [18,19], to change the hydrate pressure-temperature equilibrium conditions by injecting chemicals, such as alcohols and salts; and CO2 exchange which is most promising method to reduce the emissions of greenhouse gases [20,21].Figure 3 indicates the former three methods presented above are based on the mechanism that changes the pressure or temperature condition to destroy the stability of the original equilibrium of the hydrate deposit, resulting in its dissociation and the production of methane.However, with the relatively low energy cost and economic effectiveness, the depressurization method is generally regarded as technically feasible and the most promising one for field hydrate dissociation in the last few decades [22][23][24].In accordance with the research results published in [25], the Mallik 2002 well proved that it is feasible to use a combination of thermal stimulation and depressurization to obtain gas from permafrost [26].Recent studies have also demonstrated that using vertical wells technology to produce gas from natural hydrate reservoirs can maintain high rates for a long period in certain cases, as described by Li and Kurihara et al. [22,24].However, there is also simulation research indicating that using horizontal wells has tremendous advantages over vertical wells The following equation describes the essence of hydrate formation and dissociation and the direction of the reaction: where N H is the hydration number and N H = 5.75 is the complete hydration number, while N H = 6.0 corresponds to the average number.∆H is the enthalpy of the formation and dissociation reactions [6].
The hydration process in Equation (1) has also been thoroughly studied in the past few years as a promising method to separate methane from gaseous mixtures [8][9][10][11].Recently, several ways of developing hydrate reservoirs (Figure 2) have been put forward, including thermal stimulation [12][13][14], injection of hot brine, steam or hot water to heat the hydrate reservoir to exceed the dissociation temperature; depressurization [15][16][17], to lower the original pressure of hydrate reservoir under the dissociation pressure at a specific temperature; thermodynamic inhibitor injection [18,19], to change the hydrate pressure-temperature equilibrium conditions by injecting chemicals, such as alcohols and salts; and CO 2 exchange which is most promising method to reduce the emissions of greenhouse gases [20,21].Figure 3 indicates the former three methods presented above are based on the mechanism that changes the pressure or temperature condition to destroy the stability of the original equilibrium of the hydrate deposit, resulting in its dissociation and the production of methane.However, with the relatively low energy cost and economic effectiveness, the depressurization method is generally regarded as technically feasible and the most promising one for field hydrate dissociation in the last few decades [22][23][24].In accordance with the research results published in [25], the Mallik 2002 well proved that it is feasible to use a combination of thermal stimulation and depressurization to obtain gas from permafrost [26].Recent studies have also demonstrated that using vertical wells technology to produce gas from natural hydrate reservoirs can maintain high rates for a long period in certain cases, as described by Li and Kurihara et al. [22,24].However, there is also simulation research indicating Energies 2016, 9, 714 3 of 20 that using horizontal wells has tremendous advantages over vertical wells to obtain the gas from Class 2 and Class 3 hydrate deposits (e.g., Moridis [27] and Wang et al. [28]).However, a single well was used in all of these numerical depressurization stimulations with theoretical models.No research results about numerical depressurization stimulation with a multiple well in an actual geological CH 4 -hydrate deposit model have been reported yet.In practice, there is a tendency and need to use multi-well systems to develop a hydrate reservoir field.Therefore, it is quite interesting to simulate the dissociation process and the production behaviors in a real 3D geological CH 4 -hydrate deposit model with a multi-well system [28].
Energies 2016, 9, 714 3 of 20 to obtain the gas from Class 2 and Class 3 hydrate deposits (e.g., Moridis [27] and Wang et al. [28]).However, a single well was used in all of these numerical depressurization stimulations with theoretical models.No research results about numerical depressurization stimulation with a multiple well in an actual geological CH4-hydrate deposit model have been reported yet.In practice, there is a tendency and need to use multi-well systems to develop a hydrate reservoir field.Therefore, it is quite interesting to simulate the dissociation process and the production behaviors in a real 3D geological CH4-hydrate deposit model with a multi-well system [28].to obtain the gas from Class 2 and Class 3 hydrate deposits (e.g., Moridis [27] and Wang et al. [28]).However, a single well was used in all of these numerical depressurization stimulations with theoretical models.No research results about numerical depressurization stimulation with a multiple well in an actual geological CH4-hydrate deposit model have been reported yet.In practice, there is a tendency and need to use multi-well systems to develop a hydrate reservoir field.Therefore, it is quite interesting to simulate the dissociation process and the production behaviors in a real 3D geological CH4-hydrate deposit model with a multi-well system [28].

Status Quo in Hydrate Reservoir Modeling
Due to the complexity of the underground environment, numerical models and simulations have a great significance in the study of gas hydrate reservoir development.In recent years, various numerical simulation programs such as MH-21 HYDRES, Transport of Unsaturated Groundwater and Heat (TOUGH) + HYDRATE, HydrateResSim, STOMP-HYD and the Computer Modelling Group-Advanced Processes & Thermal Reservoir Simulator (CMG STARS) have been developed [29].
CMG-STARS which is a commercial black oil reservoir simulator can describe hydrate reservoirs.This simulator includes the Kim-Bishnoi kinetic parameters that can describe heat of dissociation and hydrate thermodynamic stability which is the core mechanism for hydrate simulation [30].Compared with other software that can be adapted to describe hydrate reservoirs, its accuracy and suitability to represent production performance from hydrate deposit in porous media has been validated by many researchers (e.g., [29,31,32]).TOUGH + HYDRATE, originated and developed at the Lawrence Berkeley National Laboratory was the first available model to simulate hydrate reservoirs, it can take in four mass components (i.e., inhibitors, water, hydrate, and methane) partitioned between four possible phases (ice, water, gas and hydrate,) and fully couple reservoir heat transfer and mass to simulate the non-isothermal dissociation process, as described by Moridis [33].The STOMP-HYD simulator uses the integral volume differencing with orthogonal grids for spatial discretization to solve the dominating conservation equations (White et al. [34]).HydrateResSim, an open-source code through the National Energy Technology Laboratories, is modified for binary and ternary hydrates, as mentioned by Garapati [35].MH-21 HYDRES was developed by the Production Method and Modeling Group to estimate production process for methane hydrate simulation, what's more, their core-test studies indicated that lowing the reservoir pressure may be the most feasible method for hydrate reservoirs at the Eastern Nankai Trough (Masuda et al. [36]).Considering the strengths and adaptability of each numerical simulation program, in this study CMG STAR was chosen as the tool to simulate the hydrate development process [37].

Objectives
The main objective of this study was to gain a better understanding of the feasible field operational factors and production processes which control the production behavior of a real 3D geological CH 4 -hydrate deposit in the Shenhu Area of the South China Sea, and to carry out a systematic sensitivity analysis of the effects of the operational parameters such as well type, well spacing, bottom hole pressure, and perforated intervals on methane recovery.Several possible production methods and combination of different operational parameters are simulated and assessed.

Geological Settings and Real Model
As shown in Figure 4, the Shenhu Area is located in the near southeast edge of the Shenhu Underwater Sandy Bench in the middle of the north slope of the South China Sea, between the Xisha Trough and the Dongsha Islands [38].Tectonically, the research area is located in the Zhu II Depression of the Pearl River Mouth Basin, which has been in a process of tectonic subsidence since the middle Miocene.In China, thegas hydrate samples were first collected at three sites (SH2, SH3, and SH7) during a recent scientific expedition conducted by the China Geological Survey in the Shenhu area of the northern South China Sea in May 2007 [39,40].Subsequently, the second Chinese marine hydrate expedition GMGS2 was performed in the Dongsha Area of the South China Sea in the summer of 2013.The hydrate resources appear to be abundant in Pearl River Mouth Basin, which has attracted much attention as a potentially important area for gas hydrate research and development in China.
logging interpretation.Using a sequential Gaussian simulation method, a porosity and permeability 3D distribution model of the hydrate-bearing layers was established.The layer porosity and hydrate saturation in the study area varies between 25%-46% and 20%-43%, respectively.As shown in Figure 5, the porosity and permeability distributions are the most important properties of the actual geological model which can reflect the core difference between the actual and ideal model about the influence of operational conditions.

Reservoir Parameters and Component Setting
This section builds a real 3D geological model on the basis of the actual geological data from the Shenhu case.Its size was 94 grids in the I-direction and 118 grids in the J-direction with a 50 m × 50 m grid, and the total area is 27.73 km 2 .In the K-direction, the model is divided unequally into 40 layers from the depth of 195 m to 229 m.Consequently, there are in total 433,680 elements.The surrounding boundaries of the model were assumed to be closed.According to the basic principle of governing equations (see Appendix A), based on currently available data from site measurements of eight wells, including water depth, thickness of the hydrate-bearing layer (HBL), the logging data, core porosity, the intrinsic permeability and seismic exploration interpretation data of the reservoir, a stochastic modeling method was employed to establish the real geological model of the target area.The top of the hydrate-bearing layers is buried 155-229 m below the seafloor (MBSF), and their thickness varies from 10 to 43 m.These hydrate layers are located at a water depth of 1108-1245 m.Porosity and permeability, which represent the physical properties of hydrate-bearing formation, were obtained from the laboratory analysis of the actual core and well logging interpretation.Using a sequential Gaussian simulation method, a porosity and permeability 3D distribution model of the hydrate-bearing layers was established.The layer porosity and hydrate saturation in the study area varies between 25%-46% and 20%-43%, respectively.As shown in Figure 5, the porosity and permeability distributions are the most important properties of the actual geological model which can reflect the core difference between the actual and ideal model about the influence of operational conditions.

Reservoir Parameters and Component Setting
This section builds a real 3D geological model on the basis of the actual geological data from the Shenhu case.Its size was 94 grids in the I-direction and 118 grids in the J-direction with a 50 m × 50 m grid, and the total area is 27.73 km 2 .In the K-direction, the model is divided unequally into 40 layers from the depth of 195 m to 229 m.Consequently, there are in total 433,680 elements.The surrounding

Reservoir Parameters and Component Setting
This section builds a real 3D geological model on the basis of the actual geological data from the Shenhu case.Its size was 94 grids in the I-direction and 118 grids in the J-direction with a 50 m × 50 m grid, and the total area is 27.73 km 2 .In the K-direction, the model is divided unequally into 40 layers from the depth of 195 m to 229 m.Consequently, there are in total 433,680 elements.The surrounding boundaries of the model were assumed to be closed.
In the simulation process of gas hydrate development, several field operational factors are considered as being relevant and necessary for the production from an actual geological CH 4 -hydrate deposit.The most important actual field operational factors which will be discussed in the subsequent sections are: well spacing (WS), bottom hole pressure (BHP), well type (WT) and perforation intervals (PI).
Table 1 summarizes the main simulation model parameters which were used in the model.Values for the actual field operational factors in the table were selected to create a reference case and subsequently discuss and analyze the spatial distribution of some parameters in the third section.Then we will use a sensitivity analysis of the actual field operational factors to discuss and put forward some personal observations about methane recovery in the fourth section.

Scenarios & Simulation
This section illustrates several scenarios about how the sensitivity analysis are made.As to the four primary actual field operational factors (well spacing, BHP, well type and perforation intervals), we establish a base case (Table 1) and then vary some factors to some degree while the other factors remain fixed at the initial (base case) value (Table 2).

Reference Case and Gas Production
The main simulation model parameters pertaining to the reference case are listed in Table 1.The combination of the field condition of reference case is 1000 m well spacing, BHP is 3000 kPa, vertical well and perforated in 4-13 layers.Figure 6 shows the evolution of the daily rate of gas production and the cumulative gas production.The 1 year production process can be divided into three periods, and the dashed lines express the borders of the different periods.Note that the daily gas production rises up rapidly to point A at around the 10th day and the rate is 753.37 ST m 3 /day, then the speed of the rise slows down and remains steady.The local maximum of gas rate daily is 4182.46ST m 3 /day before a decline begins at about point B, on the 230th day, and after that dividing point the daily gas rate tends to decline steadily.In the end, the cumulative volumetric gas production reaches about 1.091 × 10 6 ST m 3 .This tendency is attributed to a combination of: (a) a drop in the initial pressure system of the hydrate reservoir during the whole production process because of the pressure difference between the wellbore and the hydrate, and the pressure difference is the main driving force to be considered; (b) with the solid hydrate deposit decomposing and the constant driving force of depressurization, the effective permeability in the adjacent area of the well is increasing; (c) the water from the adjacent area breakthrough to the well to alleviate the effect of depressurization on gas productivity; (d) a drop in the initial temperature system resulting from the dissociation reaction would further slow the dissociation speed, especially after 230 days' production [12].
Energies 2016, 9, 714 7 of 20 well and perforated in 4-13 layers.Figure 6 shows the evolution of the daily rate of gas production and the cumulative gas production.The 1 year production process can be divided into three periods, and the dashed lines express the borders of the different periods.Note that the daily gas production rises up rapidly to point A at around the 10th day and the rate is 753.37 ST m 3 /day, then the speed of the rise slows down and remains steady.The local maximum of gas rate daily is 4182.46ST m 3 /day before a decline begins at about point B, on the 230th day, and after that dividing point the daily gas rate tends to decline steadily.In the end, the cumulative volumetric gas production reaches about 1.091 × 10 6 ST m 3 .This tendency is attributed to a combination of: (a) a drop in the initial pressure system of the hydrate reservoir during the whole production process because of the pressure difference between the wellbore and the hydrate, and the pressure difference is the main driving force to be considered; (b) with the solid hydrate deposit decomposing and the constant driving force of depressurization, the effective permeability in the adjacent area of the well is increasing; (c) the water from the adjacent area breakthrough to the well to alleviate the effect of depressurization on gas productivity; (d) a drop in the initial temperature system resulting from the dissociation reaction would further slow the dissociation speed, especially after 230 days' production [12].

Spatial Distribution of S g and S h
Figures 7 and 8 demonstrate the distribution evolution of the S g and S h in the Hydrate-Bearing Layer (HBL) during the gas production process, respectively.Figure 7 shows: (i) the waveform change of the S g distributions; (ii) the initial hydrate dissociates at around the well and expands forward to further areas; (iii) the evolution of the dissociation interface at the upper and middle region; (iv) the accumulation of gas in the upper region.Of those, (i) are results of the characteristic porosity and permeability and the heat and fluid flow; (ii) are caused by the reasonable driving force of the continuing pressure drop and spreading pressure wave; (iii) are caused by the fluids generated from the hydrate decomposition on the lower reaction front that move toward the well in the upper area; (iv) occurs because the low effective permeability k eff of the HBL in the upper area around the well [12].

Spatial Distribution of T and P
Figures and 10 demonstrate the distribution evolution of the T and P during the gas production process, respectively.Figures 9 and 10 show that: (i) once the initial temperature and pressure system is broken, the T and P around the well change more drastic and earlier than other areas; (ii) the drop of temperature and pressure in the lower region; (iii) the phenomenon of the jagged T distribution, P distribution is apparent and expanding and then turns to becomes weak with the development; (iv) the inflections come out at the interface of the dissociated and undissociated zone in the HBL.Of those, (i) and (ii) result from the velocity of pressure spreading and the characteristic of the thermal conductivity and permeability; and (iii) and (iv) result from the main advancing reaction front and the fluid flow of gas and water [12].

Spatial Distribution of T and P
Figures 9 and 10 demonstrate the distribution evolution of the T and P during the gas production process, respectively.Figures 9 and 10 show that: (i) once the initial temperature and pressure system is broken, the T and P around the well change more drastic and earlier than other areas; (ii) the drop of temperature and pressure in the lower region; (iii) the phenomenon of the jagged T distribution, P distribution is apparent and expanding and then turns to becomes weak with the development; (iv) the inflections come out at the interface of the dissociated and undissociated zone in the HBL.Of those, (i) and (ii) result from the velocity of pressure spreading and the characteristic of the thermal conductivity and permeability; and (iii) and (iv) result from the main advancing reaction front and the fluid flow of gas and water [12].Of those, (i) and (ii) result from the velocity of pressure spreading and the characteristic of the thermal conductivity and permeability; and (iii) and (iv) result from the main advancing reaction front and the fluid flow of gas and water [12].

Results and Discussion
Before discussing the different scenarios results, it is necessary to emphasize that the aim of this study was to carry out a sensitivity analysis for feasible field operational factors based on the actual 3D geological model of Shenhu area.Although the value range of the parameters may not be

Results and Discussion
Before discussing the different scenarios results, it is necessary to emphasize that the aim of this study was to carry out a sensitivity analysis for feasible field operational factors based on the actual 3D geological model of Shenhu area.Although the value range of the parameters may not be applicable to other types of hydrate reservoir, the results and discussion have contributed to unveil

Results and Discussion
Before discussing the different scenarios results, it is necessary to emphasize that the aim of this study was to carry out a sensitivity analysis for feasible field operational factors based on the actual 3D geological model of Shenhu area.Although the value range of the parameters may not be applicable to other types of hydrate reservoir, the results and discussion have contributed to unveil what the most sensitive factors are and provide guidance for actual production.Figures 11-16 present the daily gas rate, cumulative gas production and some important parameter contrasts obtained from the simulations of changing the variables one at a time.what the most sensitive factors are and provide guidance for actual production.Figures 11-16 present the daily gas rate, cumulative gas production and some important parameter contrasts obtained from the simulations of changing the variables one at a time.

Sensitivity of Producer Well Spacing
Producer spacing is an important parameter in hydrate reservoir development.As shown in Table 2, different well spacing scenarios (the well spacing ranged from 250 to 1000 m) were designed for analysis, and Figure 11 shows the cumulative gas production is extremely sensitive to the well spacing.The 1 years' production process can be divided into three periods.Before the pressure dropdown wave which is when the interactions of each well can reach an adjacent well, each scenario has an approximate trend and there is not much difference between them, which can be summarized as the initial period, Period A.
When the pressure dropdown waves of each well can interact each other at around the 10th day, the interaction with each well makes each scenario act differently until the 320th day, and we summarize this period as the interaction & storage period, which is also Period B. Although common sense would suggest that with a closer well spacing it will take the pressure dropdown wave less time to reach an adjacent well which means the pressure around the well drops faster and comes to the decomposition point easier and more decomposing gas will be produced, the results of the 250 m and 500 m well spacing scenarios show opposite trends.This phenomenon indicates that compared with pressure dropping speed, the controlling storage of each well which is reflected by the well spacing has a more powerful effect on the gas rate in this period.As for the 750 m and 1000 m scenarios, the effect of pressure interactions and controlling storage find a middle ground which makes the 750 m scenarios act best in Period B.
At approximately the 320th day, the 750 m well spacing gas rate curve and the 1000 m well spacing gas rate curves intersect with each other and after that point the 1000 m scenario produces more gas daily which indicates that the controlling storage is the conclusive factor in the long run.In the following decades the cumulative gas of 1000 m scenarios will exceed the 750 m scenario which can be predicted by the theory above and this is the so called storage period, Period C.

Sensitivity of BHP
Producer BHP is also an important parameter in the hydrate gas development.As shown in Table 2, different BHP scenarios (where the BHP ranges from 2500 to 4000 kPa) were designed for analysis.Figure 12 suggests that the cumulative gas production is extremely sensitive to the BHP, because a lower BHP results in a higher driving force, leading to a higher gas travel velocity in the HBL and in turn stimulating even more hydrate to decompose consequently.These four scenarios have the same well spacing of 1000 m, so the main force of the decomposition speed is the drawdown pressure between the hydrate reservoir pressure and the BHP, the larger the drawdown pressure is, the sooner the hydrate host formation can reach to its dissociation pressure point, and the more gas will be produced in the beginning period.Similarly, the 1 years' production time can also be divided into three periods on the basis of the daily gas production rate.The travel distance of the pressure dropdown wave is the same while the speed is quite different which is determined by the drawdown pressure between the hydrate reservoir and the BHP, so the faster the speed is, the sooner the pressure dropdown waves can reach other wells, and the sooner the interaction can happen.The initial period gas rate curves are different from those in Figure 11, and the four curves of Figure 12 are more disordered.in turn stimulating even more hydrate to decompose and the pressure to drop more in the end.For a similar reason as discussed in the former sections, the rising stage can be divided into two periods, Period A and Period B. At approximately the 260th day, the two curves intersect with each other and after that point the five vertical wells scenario produces more daily gas which indicates that the use of horizontal wells in such reservoirs may not be attractive in the long run because all layers show an upward gas migration block which is contrary to gas production from conventional gas reservoirs.

Sensitivity of Perforation Intervals of Producing Well
The perforation intervals of producing wells have an extremely significant effect on the development of hydrate gas reservoirs.Figure 15 shows that under the conditions of producer well spacing of 1000 m, bottom hole pressure of 3000 kPa, perforation intervals of 4-13 layers, the cumulative gas production and daily gas rate are the largest, and the differences between different perforation intervals are so huge that one can easily distinguish them from each other.Such huge differences may be a combined consequence by many phenomena, such as the formation pressure, non-homogeneity, temperature, grids' shape, the accuracy of program and so on.Figure 16 shows the pressure distribution of the last simulation day when perforated at different layers.It can be seen that the pressures around the perorated layers are so different that a lower pressure in the hydrate deposit produced more gas than other scenarios.The reason may be due to the fact that the upper in turn stimulating even more hydrate to decompose and the pressure to drop more in the end.For a similar reason as discussed in the former sections, the rising stage can be divided into two periods, Period A and Period B. At approximately the 260th day, the two curves intersect with each other and after that point the five vertical wells scenario produces more daily gas which indicates that the use of horizontal wells in such reservoirs may not be attractive in the long run because all layers show an upward gas migration block which is contrary to gas production from conventional gas reservoirs.

Sensitivity of Perforation Intervals of Producing Well
The perforation intervals of producing wells have an extremely significant effect on the development of hydrate gas reservoirs.Figure 15 shows that under the conditions of producer well spacing of 1000 m, bottom hole pressure of 3000 kPa, perforation intervals of 4-13 layers, the cumulative gas production and daily gas rate are the largest, and the differences between different perforation intervals are so huge that one can easily distinguish them from each other.Such huge differences may be a combined consequence by many phenomena, such as the formation pressure, non-homogeneity, temperature, grids' shape, the accuracy of program and so on.Figure 16 shows the pressure distribution of the last simulation day when perforated at different layers.It can be seen that the pressures around the perorated layers are so different that a lower pressure in the hydrate deposit produced more gas than other scenarios.The reason may be due to the fact that the upper

Orthogonal Design
With the rapid development of numerical calculations and new computational techniques in the past few decades, reservoir simulation technology has become cost-effective for model industries, not only for maximizing the production ability by factor optimization but also for assessing new production designs.It is in this study that sensitivity analysis can be conducive to assess interactions between factors, key factor determination, production forecasting, and phenomenological understanding [41].Sensitivity analysis is a tool for a model to systematically vary factors to confirm the effect of such changes [42].In practice, it has been employed in several numerical reservoir studies [41,43,44].There are two ways to put the sensitivity analysis into practice: (1) one factor at a time (OFAT); and (2) orthogonal design.OFAT is the most general means to carry out sensitivity analysis by electing a reference level for each parameter and then changing each parameter over its range while the other parameters remain unchanged [45].Although this approach can unambiguously attribute any variation that observed in the result to the single changed parameter and increase the comparability

Sensitivity of Producer Well Spacing
Producer spacing is an important parameter in hydrate reservoir development.As shown in Table 2, different well spacing scenarios (the well spacing ranged from 250 to 1000 m) were designed for analysis, and Figure 11 shows the cumulative gas production is extremely sensitive to the well spacing.The 1 years' production process can be divided into three periods.
Before the pressure dropdown wave which is when the interactions of each well can reach an adjacent well, each scenario has an approximate trend and there is not much difference between them, which can be summarized as the initial period, Period A.
When the pressure dropdown waves of each well can interact each other at around the 10th day, the interaction with each well makes each scenario act differently until the 320th day, and we summarize this period as the interaction & storage period, which is also Period B. Although common sense would suggest that with a closer well spacing it will take the pressure dropdown wave less time to reach an adjacent well which means the pressure around the well drops faster and comes to the decomposition point easier and more decomposing gas will be produced, the results of the 250 m and 500 m well spacing scenarios show opposite trends.This phenomenon indicates that compared with pressure dropping speed, the controlling storage of each well which is reflected by the well spacing has a more powerful effect on the gas rate in this period.As for the 750 m and 1000 m scenarios, the effect of pressure interactions and controlling storage find a middle ground which makes the 750 m scenarios act best in Period B.
At approximately the 320th day, the 750 m well spacing gas rate curve and the 1000 m well spacing gas rate curves intersect with each other and after that point the 1000 m scenario produces more gas daily which indicates that the controlling storage is the conclusive factor in the long run.In the following decades the cumulative gas of 1000 m scenarios will exceed the 750 m scenario which can be predicted by the theory above and this is the so called storage period, Period C.

Sensitivity of BHP
Producer BHP is also an important parameter in the hydrate gas development.As shown in Table 2, different BHP scenarios (where the BHP ranges from 2500 to 4000 kPa) were designed for analysis.Figure 12 suggests that the cumulative gas production is extremely sensitive to the BHP, because a lower BHP results in a higher driving force, leading to a higher gas travel velocity in the HBL and in turn stimulating even more hydrate to decompose consequently.These four scenarios have the same well spacing of 1000 m, so the main force of the decomposition speed is the drawdown pressure between the hydrate reservoir pressure and the BHP, the larger the drawdown pressure is, the sooner the hydrate host formation can reach to its dissociation pressure point, and the more gas will be produced in the beginning period.
Similarly, the 1 years' production time can also be divided into three periods on the basis of the daily gas production rate.The travel distance of the pressure dropdown wave is the same while the speed is quite different which is determined by the drawdown pressure between the hydrate reservoir and the BHP, so the faster the speed is, the sooner the pressure dropdown waves can reach other wells, and the sooner the interaction can happen.The initial period gas rate curves are different from those in Figure 11, and the four curves of Figure 12 are more disordered.
The B period can be called the speeding period because the differences of the curves are so obvious and all the gas rates are increasing which indicates the drawdown pressure effect in this period is the dominant factor and the 2500 kPa curve areis rising rapidly.
At approximately the 250th day, all the gas rate curves are declining and have a trend to coincide with each other when the gap between the hydrate reservoir pressure and the BHP becomes smaller and even vanishes.In the following decades, when the gap is gone, the reaction will come to a phase equilibrium again and no more extra gas will be produced.The accumulated gas of the 2500 kPa scenarios will no doubt be the most because the initial gap is the largest which can be predicted by the theory above and this is the so called declining period, Period C.

Sensitivity of Well Types
Well type is also an important parameter in hydrate gas development.As shown in Table 2, different well type scenarios were designed for the analysis, and Figure 13 shows that the distance between a horizontal well's head and end is equal to the well spacing of vertical wells.Figure 14 indicates that the daily gas production is extremely sensitive to the well type at the beginning, because a horizontal well has a more effective contact area, resulting in more controlled storage and in turn stimulating even more hydrate to decompose and the pressure to drop more in the end.For a similar reason as discussed in the former sections, the rising stage can be divided into two periods, Period A and Period B. At approximately the 260th day, the two curves intersect with each other and after that point the five vertical wells scenario produces more daily gas which indicates that the use of horizontal wells in such reservoirs may not be attractive in the long run because all layers show an upward gas migration block which is contrary to gas production from conventional gas reservoirs.

Sensitivity of Perforation Intervals of Producing Well
The perforation intervals of producing wells have an extremely significant effect on the development of hydrate gas reservoirs.Figure 15 shows that under the conditions of producer well spacing of 1000 m, bottom hole pressure of 3000 kPa, perforation intervals of 4-13 layers, the cumulative gas production and daily gas rate are the largest, and the differences between different perforation intervals are so huge that one can easily distinguish them from each other.Such huge differences may be a combined consequence by many phenomena, such as the formation pressure, non-homogeneity, temperature, grids' shape, the accuracy of program and so on.Figure 16 shows the pressure distribution of the last simulation day when perforated at different layers.It can be seen that the pressures around the perorated layers are so different that a lower pressure in the hydrate deposit produced more gas than other scenarios.The reason may be due to the fact that the upper layer has a lower initial pressure and it is easy and fast to decline to its decomposition pressure point with constant BHP.Therefore, perforation intervals have an extremely strong effect on the production.The production well perforation interval should be located in the upper deposit in the development of the Shenhu case.

Orthogonal Design
With the rapid development of numerical calculations and new computational techniques in the past few decades, reservoir simulation technology has become cost-effective for model industries, not only for maximizing the production ability by factor optimization but also for assessing new production designs.It is in this study that sensitivity analysis can be conducive to assess interactions between factors, key factor determination, production forecasting, and phenomenological understanding [41].Sensitivity analysis is a tool for a model to systematically vary factors to confirm the effect of such changes [42].In practice, it has been employed in several numerical reservoir studies [41,43,44].There are two ways to put the sensitivity analysis into practice: (1) one factor at a time (OFAT); and (2) orthogonal design.OFAT is the most general means to carry out sensitivity analysis by electing a reference level for each parameter and then changing each parameter over its range while the other parameters remain unchanged [45].Although this approach can unambiguously attribute any variation that observed in the result to the single changed parameter and increase the comparability of the results, it cannot take the interactions between the parameters which are due to the failure of one parameter to have the same effect on the response while another parameter is changed too [42].
In this paper, although the analysis of field operational conditions has been made by the OFAT in the last few sections, we still cannot draw a specific conclusion about which is the most sensitive parameter and the effect order among well spacing, bottom hole pressure and perforation intervals.Thus, we employ an orthogonal design to determine the relative importance of field operation variables on the cumulative gas produced from an actual geological CH 4 -hydrate deposit.Figure 17 indicates the production performance prediction of the orthogonal design experiment.
Energies 2016, 9, 714 14 of 20 of the results, it cannot take the interactions between the parameters which are due to the failure of one parameter to have the same effect on the response while another parameter is changed too [42].
In this paper, although the analysis of field operational conditions has been made by the OFAT in the last few sections, we still cannot draw a specific conclusion about which is the most sensitive parameter and the effect order among well spacing, bottom hole pressure and perforation intervals.Thus, we employ an orthogonal design to determine the relative importance of field operation variables on the cumulative gas produced from an actual geological CH4-hydrate deposit.Figure 17 indicates the production performance prediction of the orthogonal design experiment.
Table 3 shows the basic data of the orthogonal design experiment and the analysis results.Mean A, B and C are taken the average value of cumulative gas in three groups of experiment, and range is the difference between the maximum and the minimum of the mean for a certain factor.The effect of factors can be analyzed by comparing the numerical value of a range, and the larger the range, the greater effect.The one factor at a time design and the orthogonal design indicate that the order of the effects of the factors on gas yield was perforation intervals > bottom hole pressure > well spacing, and the best combination of operation condition is A2, B1and C1 which is well spacing 750 m, BHP 2500 kPa and perforation in the 4-13 layers, whereby the cumulative gas production can be 1,309,742 m 3 .Table 3 shows the basic data of the orthogonal design experiment and the analysis results.Mean A, B and C are taken the average value of cumulative gas in three groups of experiment, and range is the difference between the maximum and the minimum of the mean for a certain factor.The effect of factors can be analyzed by comparing the numerical value of a range, and the larger the range, the greater effect.The one factor at a time design and the orthogonal design indicate that the order of the effects of the factors on gas yield was perforation intervals > bottom hole pressure > well spacing, and the best combination of operation condition is A2, B1and C1 which is well spacing 750 m, BHP 2500 kPa and perforation in the 4-13 layers, whereby the cumulative gas production can be 1,309,742 m 3 .

Conclusions
The following conclusions can be taken from this study: (1) As for the actual geological model of the Shenhu area, the gas accumulation of the basis scenario of operational conditions with 1000 m well spacing, 3000 kPa BHP and perforated 4-13 intervals can reach 1.02 × 10 6 S Tm 3 after one years' production, and the peak daily gas rate can reach approximately 4000 STm 3 /day.That outcome can be used as reference data when it comes to the real field development.
(2) The operation conditions indeed have an important influence on the production performance.
The well spacing analysis indicates that the total production period can be divided into three periods according to the degrees of the influencing factors that are initial period, interacting & storage period and storage period.Different well spacing acts similarly in period A because of the pressure interaction has not happened, the pressure dropdown waves and the controlling storage work together in the period B, and the controlling storage plays the main role in the following period when the pressure interaction is slight after the long production time.
(3) The BHP analysis indicates that the total production period can also be divided into three periods according to the size of pressure gap that are initial period, speeding period and declining period.Different BHP reflects different pressure gaps between the hydrate reservoir and the BHP, that is the driving force which subsequently results in a higher gas travel velocity in the host formation which subsequently in turn stimulates even more hydrate dissociation.In the following decades, when the gap is gone, the reaction will remain in phase equilibrium and no more extra gas will be produced.(4) The perforation intervals analysis indicates that the pressure around the perforated layers is so dominant that lower pressure in the hydrate deposit produced more gas than other scenarios.The perforation intervals have an extremely strong effect on the production.The production well perforation interval should be located in the upper deposit in the development of the Shenhu case.(5) The well type analysis indicates the use of horizontal wells in such reservoirs may not be attractive in the long run as all layers showed an upward gas migration block which is contrary to gas production from conventional gas reservoirs.(6) One factor at a time design and orthogonal design indicate that the order of the effects of the factors on gas yield was perforation intervals > bottom hole pressure > well spacing.
ρ g gas density (g/cm 3 ) ρ w water density (g/cm 3 ) ρ H hydrate density (g/cm 3 ) v gx gas velocity along the x-axis (m/s) v gy gas velocity along the y-axis (m/s) v gz gas velocity along the z-axis (m/s)

Figure 2 .
Figure 2. The mechanism of the three main methods to develop hydrate reservoirs.(A) Thermal Stimulation; (B) Depressurization; (C) Inhibitor Injection.

Figure 3 .
Figure 3.The mechanism for changing the pressure or the temperature conditions to destroy the stability of the original equilibrium and cause a phase change.

Figure 2 .
Figure 2. The mechanism of the three main methods to develop hydrate reservoirs.(A) Thermal Stimulation; (B) Depressurization; (C) Inhibitor Injection.

Figure 2 .
Figure 2. The mechanism of the three main methods to develop hydrate reservoirs.(A) Thermal Stimulation; (B) Depressurization; (C) Inhibitor Injection.

Figure 3 .
Figure 3.The mechanism for changing the pressure or the temperature conditions to destroy the stability of the original equilibrium and cause a phase change.Figure 3. The mechanism for changing the pressure or the temperature conditions to destroy the stability of the original equilibrium and cause a phase change.

Figure 3 .
Figure 3.The mechanism for changing the pressure or the temperature conditions to destroy the stability of the original equilibrium and cause a phase change.Figure 3. The mechanism for changing the pressure or the temperature conditions to destroy the stability of the original equilibrium and cause a phase change.

Figure 4 .
Figure 4.A map showing the location of the study area, drilling sites, and the confirmed gas hydrate distribution in the Shenhu area.

Figure 5 .
Figure 5.The distribution of permeability of an actual geological model on the northern continental slope of the South China Sea.

Figure 4 .
Figure 4.A map showing the location of the study area, drilling sites, and the confirmed gas hydrate distribution in the Shenhu area.
the reservoir, a stochastic modeling method was employed to establish the real geological model of the target area.The top of the hydrate-bearing layers is buried 155-229 m below the seafloor (MBSF), and their thickness varies from 10 to 43 m.These hydrate layers are located at a water depth of 1108-1245 m.Porosity and permeability, which represent the physical properties of hydrate-bearing formation, were obtained from the laboratory analysis of the actual core and well logging interpretation.Using a sequential Gaussian simulation method, a porosity and permeability 3D distribution model of the hydrate-bearing layers was established.The layer porosity and hydrate saturation in the study area varies between 25%-46% and 20%-43%, respectively.As shown in Figure5, the porosity and permeability distributions are the most important properties of the actual geological model which can reflect the core difference between the actual and ideal model about the influence of operational conditions.

Figure 4 .
Figure 4.A map showing the location of the study area, drilling sites, and the confirmed gas hydrate distribution in the Shenhu area.

Figure 5 .
Figure 5.The distribution of permeability of an actual geological model on the northern continental slope of the South China Sea.

Figure 5 .
Figure 5.The distribution of permeability of an actual geological model on the northern continental slope of the South China Sea.

Figure 6 .
Figure 6.The volumetric rate of daily gas production and the cumulative gas production of reference case in one years' simulation.

Figures 7
Figures 7 and 8 demonstrate the distribution evolution of the Sg and Sh in the Hydrate-Bearing Layer (HBL) during the gas production process, respectively.Figure 7 shows: (i) the waveform change of the Sg distributions; (ii) the initial hydrate dissociates at around the well and expands forward to further areas; (iii) the evolution of the dissociation interface at the upper and middle region; (iv) the accumulation of gas in the upper region.Of those, (i) are results of the characteristic porosity and permeability and the heat and fluid flow; (ii) are caused by the reasonable driving force of the continuing pressure drop and spreading pressure wave; (iii) are caused by the fluids generated from the hydrate decomposition on the lower reaction front that move toward the well in the upper area; (iv) occurs because the low effective permeability keff of the HBL in the upper area around the well [12].
Figures 7 and 8 demonstrate the distribution evolution of the Sg and Sh in the Hydrate-Bearing Layer (HBL) during the gas production process, respectively.Figure 7 shows: (i) the waveform change of the Sg distributions; (ii) the initial hydrate dissociates at around the well and expands forward to further areas; (iii) the evolution of the dissociation interface at the upper and middle region; (iv) the accumulation of gas in the upper region.Of those, (i) are results of the characteristic porosity and permeability and the heat and fluid flow; (ii) are caused by the reasonable driving force of the continuing pressure drop and spreading pressure wave; (iii) are caused by the fluids generated from the hydrate decomposition on the lower reaction front that move toward the well in the upper area; (iv) occurs because the low effective permeability keff of the HBL in the upper area around the well [12].

Figure 6 .
Figure 6.The volumetric rate of daily gas production and the cumulative gas production of reference case in one years' simulation.

Figure 7 .Figure 8 .
Figure 7.The evolution of spatial distribution of Sg during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

Figure 7 .Figure 7 .Figure 8 .
Figure 7.The evolution of spatial distribution of S g during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

Figure 8 .
Figure 8.The evolution of spatial distribution of S h during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

3. 3 .
Figures 9 and 10 demonstrate the distribution evolution of the T and P during the gas production process, respectively.Figures9 and 10show that: (i) once the initial temperature and pressure system is broken, the T and P around the well change more drastic and earlier than other areas; (ii) the drop of temperature and pressure in the lower region; (iii) the phenomenon of the jagged T distribution, P distribution is apparent and expanding and then turns to becomes weak with the development; (iv) the inflections come out at the interface of the dissociated and undissociated zone in the HBL.

Figure 9 .Figure 10 .
Figure 9.The evolution of spatial distribution of T during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

Figure 9 .Figure 9 .Figure 10 .
Figure 9.The evolution of spatial distribution of T during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

Figure 10 .
Figure10.The evolution of spatial distribution of P during gas production from the methane hydrate deposit at some site of the Shenhu Area, South China Sea.

Figure 11 .
Figure 11.The volumetric rate of daily gas production and the cumulative gas production of four scenarios with different well spacing in one years' simulation.

Figure 11 .
Figure 11.The volumetric rate of daily gas production and the cumulative gas production of four scenarios with different well spacing in one years' simulation.

Figure 12 .
Figure 12.The volumetric rate of daily gas production and the cumulative gas production of four scenarios with different BHP in one years' simulation.

Figure 12 .
Figure 12.The volumetric rate of daily gas production and the cumulative gas production of four scenarios with different BHP in one years' simulation.

Figure 13 .
Figure 13.The pressure contrast of final moments with different well type.

Figure 14 .
Figure 14.The volumetric rate of daily gas production and the cumulative gas production of two scenarios with different well type in one years' simulation.

Figure 13 .
Figure 13.The pressure contrast of final moments with different well type.

Figure 13 .
Figure 13.The pressure contrast of final moments with different well type.

Figure 14 .
Figure 14.The volumetric rate of daily gas production and the cumulative gas production of two scenarios with different well type in one years' simulation.

Figure 14 .
Figure 14.The volumetric rate of daily gas production and the cumulative gas production of two scenarios with different well type in one years' simulation.

Figure 15 .Figure 15 .
Figure 15.The volumetric rate of daily gas production and the cumulative gas production of two scenarios with different perforation intervals in one years' simulation.

Figure 15 .Figure 16 .
Figure 15.The volumetric rate of daily gas production and the cumulative gas production of two scenarios with different perforation intervals in one years' simulation.

Figure 16 .
Figure 16.The pressure contrast of final moments with different perforation intervals.

Figure 17 .
Figure 17.The cumulative gas production of the orthogonal design experiment.

Figure 17 .
Figure 17.The cumulative gas production of the orthogonal design experiment.

Table 1 .
Summary of main simulation model parameters were used in the model.

Table 2 .
Summary of the changing factor and its range.

Table 3 .
The basic data of the orthogonal design experiment and sensitive factor analysis.L (3 3 ).

Table 3 .
The basic data of the orthogonal design experiment and sensitive factor analysis.L (3 3 ).