Laboratory Study on Hydrate Production Using a Slow, Multistage Depressurization Strategy

Key Laboratory of Gas Hydrate, Ministry of Natural Resources, Qingdao Institute of Marine Geology, Qingdao 266071, China Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266237, China College of Petroleum Engineering, China University of Petroleum, Qingdao 266555, China Petroleum Systems Engineering, Faculty of Engineering and Applied Science, University of Regina, Regina, Canada


Introduction
Natural gas hydrate (NGH) is an ice-like compound formed from water and gas molecules under relatively low-temperature and high-pressure conditions [1] NGH is considered as the most promising alternative fossil fuel due to its great energy potential [2]. Hydrate-related field exploration and exploitation attract attentions from both the international ocean discovery program (IODP) and national plans such as the United States, Japan, China, South Korea, and India [3,4].
Depressurization is also restricted by potential geohazards such as severe sand production [14][15][16], wellbore instability [17], and deformation of the sediment [18,19]. These potential geohazards are mainly controlled by the mechanical properties of HBS [20,21]. Both numerical simulation [22] and experiments [23] reveal that sand production rate tends to become severe with the increase in the magnitude of total pressure-drop and depressurizing rate. Hydrate production test can be forced to terminate once the sand production rate exceeds the tolerance of wellbore or platform [10]. Besides, the influence of geological parameters, such as compressibility [24,25], absolute permeability [26,27], relative permeability [25,28], and clayey content [29], on the gas production behavior has been covered both in lab-scale [7,9] and field-scale [10] studies.
The depressurization parameters, which mainly include pressure-drop magnitude and depressurizing rate, would also affect gas production behaviors [9,[30][31][32]. Aggressive depressurization would increase gas production rate in a certain extent [5,33], but causes unpredictable occurrence of hydrate reformation [34]. There are many pathways to reach the same pressure-drop magnitude such as one-step quick depressurization and slow-stepwise depressurization [30]. The number of steps, the pressure-drop magnitude in each step, and the average depressurizing rate could be used to characterize the depressurization pathways under the same total pressure-drop magnitude.
Therefore, balancing the gas productivity and potential geohazards is crucial during depressurization [35,36]. The influences of depressurization pathways on gas-water production behaviors and formation response need to be clarified. In this paper, we form methane hydrate in a flexible rubber sleeve. Two depressurizing gas production experiments were conducted at constant confining pressure and temperature. A numerical simulation model based on CMG-STARS module was developed to match the measured gas production, as well as pressure evolutionary behaviors. Then, the verified numerical model was employed to investigate the influence of depressurization pathways on gas-water production behaviors, as well as formation pressure-temperature (P-T) and hydrate saturation responses.

Experiments
The hydrate production simulation device is shown in Figure 1. The device consists of (a) a self-developed cylindri-cal pressure vessel system, (b) a reactant supply system for injecting gas and water into the pressure vessel, (c) a thermostatic chamber to monitor the temperature, and (d) a system controlling and data acquisition system.
The main body of the pressure vessel is made from titanium alloy with maximum experimental pressure of 15 MPa. A flexible rubber sleeve (Φ60 mm × 1000 mm) is installed inside the titanium alloy vessel. The annulus between the rubber sleeve and inner of the titanium alloy vessel is filled with confining oil. Sediment for HBS formation is packed inside the rubber sleeve. Five pressure sensors (P 1~P5 ), with measurement accuracy of 0.01 MPa and measuring ranges of 10 MPa, are installed linearly on the rubber sleeve to test the real-time pore pressure. The distance from P 1 to the outlet of the vessel is 50 mm, and the distance between each of the two sensors is 200 mm.
The maximum effective confining pressure (difference between confining pressure and pore pressure) for the rubber sleeve is 8 MPa when it was packed with sediment. During hydrate formation and dissociation, the temperature and pressure of the confining oil remain constant, indicating that a constant temperature and constant pressure boundary condition could be reached during the experiments. Depressurization process is achieved by adjusting the back pressure valves at the outlet.
Quartz sand with diameter of 0.12 mm~0.18 mm was used in the experiments. The average porosity of the sediment is 39.9%. Methane gas with purity of 99.99% and distilled water were used to form hydrate. Detailed sand packing and fluid injection procedures could be found in Jin et al. [37].
We tested the intrinsic permeability of the sediment before hydrate formation. The confining pressure was adjusted to 2 MPa after the sediment was packed into the rubber sleeve. Then, distilled water was injected into the sediment at a rate of 35 ml/min. During the process of water injection, the pressure in the sediment was recorded by a pressure sensor. When the pressure values remained constant, we considered that the water flow in the sediment is a steady-state flow. The pressure difference between P 1 and P 5 was about 0.11 MPa, providing that the viscosity of water at 300 K is 0.855 mPa·s. The intrinsic permeability of the Syringe pump Back pressure valve Back pressure controller Figure 1: Schematic diagram of the hydrate production simulation device.
2 Geofluids sediment could be obtained according to Darcy's law, which is about 606 mD in this manuscript. The average hydrate saturation in the HBS, which is determined by the volume of gas consumption, was controlled around 30%. The confining pressure was set as 8 MPa, and the temperature was set as 3.5°C while forming HBS. Once the HBS formation process was completed, the confining pressure was adjusted to 6 MPa, while the effective confining pressure was about 2 MPa.
The one-step quick depressurization and two-step multistage depressurization schedules were applied to simulate hydrate production process. In the one-step quick depressurization experiment, the pressure at the outlet was controlled to decrease to the atmospheric pressure within 5 minutes. The pressure-drop magnitude was 4 MPa during gas production. In the two-step depressurization experiment, the pressure at the outlet was firstly controlled at 2 MPa for 4 hours. Then, it was decreased to the atmospheric pressure within 5 minutes.

Numerical Modeling and Verification
The TOUGH+HYDRATE (Lawrence Berkeley National Laboratory, USA) [38], MH21-Hydrates (AIST Group, Japan) [39], STOMP-HYD (Pacific Northwest National Laboratory, USA) [39], and CMG-STARS (Computer Modeling Group Ltd., Canada) [39] are widely adopted to simulate hydrate production processes. In these simulators, the hydrate intrinsic dissociation constant is normally taken from bulk hydrate dissociation. However, recent study documented that the hydrate intrinsic dissociation constant varies with the characteristics of porous media [40]. Therefore, numerical simulation should firstly verify the hydrate intrinsic dissociation constant based on experimental results.
The CMG-STARS is primarily developed for thermal recovery. The simulator considers the coupling effect of fluid flow, heat transfer, and fines transport in porous media and wellbore. In simulating hydrate production process, the hydrate is considered as "unmovable heavy oil" phase filled in the porous media [40,41], while the other part of the pore space is saturated by either water or gas phase. Based on Kim-Bishnoi model, Uddin et al. [42] modified the hydrate dissociation kinetic equation into the following equation: where dc h /dt is the hydrate dissociation rate in gmol/ [m 3 ·d]. k 0 d represents the intrinsic hydrate decay rate in gmol/[m 2 ·kPa·d]. A h is the specific area of hydrate particles, 3 × 10 5 m 2 /m 3 . S w and S h are water saturation and hydrate saturation, respectively. ϕ is the porosity of the sediment, ϕ = 39:9%. ΔE is the activation energy, 76.516 kJ/[mol·K]. R is the universal gas constant, 8.3144 J/[mol·K]. T is the temperature in K. P eq is the water-hydrate-vapour equilibrium pressure in kPa. KðP, TÞ = P e /P g .
In Equation (1), the mass concentration of water and hydrate during dissociation can be written as D w = ϕS w ρ w , D h = ϕS h ρ h . Equation (1) can be written as follows: where K D is the modified hydrate intrinsic dissociation constant for CMG-STARS.
To simulate the constant temperature and constant pressure boundary experimental conditions, a numerical model is built with three-dimensional cylindrical shape ( Figure 2). The cross-sectional area and length of the model are the same with the rubber sleeve (Φ60 mm × 1000 mm). The model is divided into 20 grid blocks in the I-direction and 25 grid blocks in the R-direction. The cross section of the model is divided into 20 grid blocks in the θ -direction. Each grid block has the same lengths in the I-direction (50 mm) and R-direction (3 mm), and the same coverage angle of 18 degree in the θ -direction. In order to investigate physical response inside the formation, a longitudinal section (a-a in Figure 3) was taken along the I-direction. The longitudinal section crosscuts the model in the R-direction. There is no mass exchange between the confining oil and the HBS.

Geofluids
The initial average hydrate saturation is 30% according to the experiments. The input parameters for numerical simulation are listed in Table 1. Detailed parameters of adjusting and regressing methods can be found in Ajayi et al. [41]. The one-step depressurization experiment is used for matching the modified hydrate intrinsic dissociation constant. The two-step depressurization experiment is used for verifying the adaptability of the numerical model.
Relative permeability of the HBS was calculated from Stone's equation [43].
where k rA and k rG are the relative permeabilities of water and gas, respectively. S irA is the immobile water saturation, whereas S irG is the immobile gas saturation. S irA = 0:25, S irG = 0:01, n = 3:5, and nG = 2:5 [44]. Comparisons between the experimental and numerical results for one-step quick depressurization and two-step multistage depressurization are shown in Figures 3 and 4, respectively. Solid lines in Figures 3 and 4 represent the numerical results, while the dash lines represent the experimental results. Figures 3 and 4 indicate that the gas production behaviors between numerical and experimental results match favorable with each other. The one-step quick depressurization witnessed only one local peak in transient gas production rate during the whole production process. However, the two-step multistage pressurization procedure saw another local peak in transient gas production rate immediately after the second depressurization step.
Under the same initial average hydrate saturation, the same boundary temperature, the same pore pressure, and the same host sediment, the cumulative gas volume from two-step multistage depressurization would be a little more than that from one-step quick depressurization. However,   the duration for gas production process from two-stepwise depressurization is shorter than that from one-step quick depressurization.
The pressure evolutionary behaviors from the experiments showed obvious hysteresis, though the trend is similar with that from numerical simulation. This could be attributed to the reformation of hydrate or icing at the inlet of pressure sensor in the experiment [45]. Once the hydrate particles accumulate at the inlet of the pressure sensor, the pressure sensor would not be able to record the real-time pore pressure. Further depressurization would dissociate hydrate accumulated at the inlet of pressure sensor; hence, a plummet in pressure would follow. The modified hydrate intrinsic dissociation constant K D is 2:43 × 10 4 mol/m 2 · Pa according to the regression.

Depressurization Pathways.
Based on the verified numerical model, field hydrate reservoir temperature and pressure conditions in 2017 offshore methane hydrate production test in the Nankai Trough of Japan are taken as the initial conditions. The average pore pressure and temperature at the middle of the hydrate reservoir are 13 MPa and 13.5°C, respectively [33]. The confining pressure of the hydrate reservoir is set as 15 MPa according to hydrostatic pressure, and the confining temperature is 17°C. Three hydrate saturation conditions (~75%,~50%, and~10%) would be considered. The size and boundary conditions of the model are the same with the experiments (Table 1).
During depressurization, the driving force for hydrate dissociation comes from the difference between waterhydrate-vapour equilibrium pressure and real-time bottom-hole flow pressure. Therefore, the water-hydrate-vapour equilibrium pressure at the reservoir temperature (10.0 MPa at 13.5°C) is taken as the reference to choose the depressurization pathways. The total pressure drop is fixed to 4 MPa, which is the same with the experiments. Hence, the final pore pressure for all simulation cases is 6 MPa. We defined six depressurization pathways to reach the total pressure-drop magnitude of 4MPa ( Figure 5).
In Figure 5, the 1 st pathway represents one-step quick depressurization, while paths 2~6 are multistage depressurization pathways. In the paths 2 nd~4th , the bottom-hole pressure is decreased to 6 MPa within 3 hours, indicating the same average depressurizing rate of 4 MPa/3 h. The  5 Geofluids depressurization processes are divided into 2 steps, 4 steps, and 16 steps equidistantly in the paths 2 nd~4th , respectively. In the 5 th path, the bottom-hole pressure is decreased to 6 MPa within 4 hours, indicating an average depressurizing rate of 4 MPa/4 h. The 6 th path is taken as an example to investigate the influence of shut-in operation during production. The depressurization procedure for the 6 th path is the same with the 2 nd pathway except the shut-in operation (from 1 h to 2 h).

Cumulative
Gas-Water Production. Figure 6 shows the evolutionary behaviors of cumulative gas and water when the initial hydrate saturation is 10% (the initial water saturation is 90%). It could be concluded from Figure 6 that depressurization pathways could dramatically affect the final cumulative gas and water volume, even though the total pressure-drop magnitude is fixed. Compared with the onestep quick depressurization (the 1 st path), stepwise depressurization could increase the final cumulative gas volume and decrease the final cumulative water volume to some extent. This implies that depressurization mode would have some significant influences on formation multiphase seepage field, as well as the residual fluid saturation. Therefore, we define the ratio of final cumulative gas volume and water volume as the integrated gas-to-water ratio (IGWR, Equation (4)). The IGWR could be used to evaluate the influence of depressurization pathways on fluid production behaviors.
When the total pressure-drop magnitude is fixed, we use the depressurizing magnitude in each step and the number of steps to identify and characterize different depressurization pathways. Under the same total pressure drop, the influence of depressurizing magnitude in each step and the number of steps on IGWR are shown in Figure 7. Generally, the IGWR increases greatly with the increase in initial hydrate saturation. Under the same initial hydrate saturation, the IGWR decreases linearly with the increase in depressurizing magnitude in each step, while increases logarithmically with the increase in the number of steps.
In Figure 7(a), if the linear fitting curves for the IGWR and depressurizing magnitude in each step are extended leftwards to the y-axis, we obtain the theoretical maximum IGWRs for stepwise depressurization. The maximum IGWR is mainly affected by the initial hydrate saturation (Figure 8). The theoretical maximum IGWR could be reached only if the depressurizing magnitude in each step becomes minimum, indicating that the depressurization procedure is conducted linearly at a constant depressurizing rate. However, this is almost impossible to be realized in a field operation.
On the other hand, when the number of steps is low, the IGWR increases rapidly with the increase in the number of steps. However, the IGWR keeps almost constant once the number of steps exceeds a certain value. The logarithmic relationships between IGWR and the number of steps indicate that it is unnecessary to maximize the number of steps in fieldwork. In short, stepwise depressurization is important for increasing the IGWR at the same total magnitude of pressure drop. However, it is unnecessary to maximize the number of steps, which would affect the difficulties of field operation.
A comparison between the 5 th path and the 4 th path indicates that a slower average depressurizing rate would improve the IGWR within a certain scope, although the depressurizing magnitude in each step and the number of steps are the same. However, the IGWR for the 2 nd path almost equals to that for the 6 th path, indicating that the  6 Geofluids shut-in operation could barely affect the IGWR. When the average depressurizing rate is low, the transient hydrate dissociation rate would be relatively low. Hence, gas phase would be separated from the water phase. Under the condition of relatively low seepage rate, the water phase cannot be carried by gas phase. As a result, the IGWR increases with the decrease in the average depressurizing rate.

Temperature Evolutionary
Behavior. The temperature of the HBS during hydrate production is mainly determined by the combination of endothermic hydrate dissociation process, Joule-Thomson effect, and heat supply from the confining formation. The 3 rd path, with an initial hydrate saturation of 75%, is taken as an example to explore temperature field changing behaviors. Figure 9 shows the temperature distribution on the section a-a (Figure 3). Point A in Figure 9 represents the geometric core of the section a-a. The right side of the model (marked by rightwards arrow) represents the outlet of the model. The low-temperature area initiates at the outlet (production well) and propagates in a nonpiston style towards the far  7 Geofluids side of the model in the I-direction. However, the temperature field shows obvious axisymmetrical evolutionary behaviors in the R-direction once the hydrate dissociation process is initiated. We attribute the axisymmetrical characteristics to the cylindrical heat supply from the confining oil. These evolutionary phenomena in temperature field were also observed in other simulation cases.
To further investigate the influence of depressurization mode on temperature, the depressurization pathways in Figure 5 are divided into two categories. The average depressurizing rate for the 2 nd~4th paths is 4 MPa/3 h, while the numbers of steps and pressure drop in each step are different. Therefore, the 2 nd~4th paths can be used to determine the influence of depressurization mode on temperature under the condition of the same total pressure-drop magnitude and the same average depressurizing rate. On the other hand, the average depressurizing rates for the 1 st , 4 th , and 5 th paths are 4 MPa/5 min, 4 MPa/3 h, and 4 MPa/4 h, respectively. Hence, they could be set into the same category to explore the influence of overall depressurizing rate.
Temperature of the confining oil is constant and equivalent for all simulation cases. Therefore, temperature differences among different simulation cases are probably caused by the combination of hydrate dissociation and Joule-Thomson effect, which are both endothermic processes. Figure 10 shows that the minimal temperature value decreases with the increase in hydrate saturation. It is easy to understand because higher hydrate saturation implies more heat would be consumed by hydrate dissociation.
Most interestingly, it can be seen from Figures 10(a), 10(c), and 10(e) that a decrease in the overall depressurizing rate could significantly alleviate the temperature drop. This is in consistency with that obtained by Zhao et al. [46]. Furthermore, a slower depressurization mode would postpone the appearance of minimal temperature value to a certain extent.
Under the same total pressure-drop value (4 MPa) and average depressurizing rate (4 MPa/3 h), increase in the number of steps (or decrease in the pressure-drop magnitude in each step) would benefit for the smoothness of the temperature behaviors. The minimal temperature value decreases in a very narrow range with the increase in the number of steps. However, the depressurization modes have little influences on temperature rebounding process once the production duration exceeds the minimal temperature.
In field operation, temperature drop is an unwanted factor during hydrate production because of potential hydrate reformation risk [34]. However, we infer from the results that we could hardly avoid temperature drop by only increasing the number of steps or decreasing the magnitude of pressure drop in each step. The best choice for avoiding temperature drop would be decreasing the average depressurizing rate. A slow, multistage depressurization strategy would be the best choice.

Hydrate Saturation Evolutionary
Behavior. The hydrate saturation field evolutionary behaviors for the 3 rd depressurization pathway when initial hydrate saturation is 75% are shown in Figure 11. Figure 11 suggests that the hydrate dissociation initiates at the edge of the model near the outlet (t = 0:0 h) and then propagates uniformly and randomly (t = 0:5 h). This would be probably caused by the abundant heat supply from the confining oil.
Hereafter, the hydrate dissociation front extends towards the far side of the model in the I-direction. The hydrate dissociation front showed obvious nonpiston characteristics, which corresponds to the temperature field changing behaviors in Figure 9. In the R-direction, the hydrate saturation field shows obvious axisymmetric evolution behaviors, indicating that the hydrate dissociation process is mainly controlled by the temperature field.
To determine the influence of depressurization modes on hydrate dissociation processes, we define the hydrate dissociation duration as the time when the hydrate saturation at the core of the far side of model (point B in Figure 11) is decreased to less than 0.1%. We collected the hydrate dissociation duration for all simulation cases under different hydrate saturations. The results are shown in Figure 12.
It could be concluded from Figure 12 that hydrate dissociation duration varies greatly with depressurization mode. The hydrate dissociation duration would increase remarkably with the decrease of average depressurizing rate (paths   (Figure 2) during the 3 rd depressurization when initial hydrate saturation is 75%. Point A represents the geometric core of the section a-a. 8 Geofluids 1, 4, and 5). However, when the average depressurizing rate is fixed, the depressurization mode could hardly affect the hydrate dissociation duration (paths 2, 3, and 4). However, the shut-in operation (path 6) would postpone the hydrate dissociation duration greatly (compared with path 2).

4.5.
Influence of Well Shut-In Operation. During marine hydrate production process, the gas production process might have to be suspended due to failure in down-hole equipment or irresistible factors such as typhoon [13]. Therefore, shut-in operation is unavoidable in fieldwork, especially for a long-term production. The temperature and pressure 52% 40% 29% 18% Figure 11: The hydrate evolutionary behaviors during the 3rd depressurization path when initial hydrate saturation is 75%. Point B represents the center of the far side of mode on the section a-a.    Figure 13: Formation pore pressure-temperature response to the shut-in operation. 10 Geofluids recovery behaviors at point A ( Figure 9) are taken to explore the influence of shut-in operation ( Figure 13). It could be concluded from Figure 13 that the pore pressure would rebound logarithmically once the gas production is suspended. The higher the hydrate saturation, the faster the pressure recovery rate would be. However, at the late stage of shut-in period, the pressure for the lower hydrate saturation (10%) case exceeds the higher hydrate saturation (50% or 75%) cases. This is because the temperature of the sediment for lower hydrate saturation (10%) case is much higher than that when hydrate saturation is 50% and 75%.
When the shut-in operation is ended, the pore pressure would decrease to the regular level (path 2) instantaneously, whereas the temperature decreases logarithmically. This is quite important for field operation because it can be inferred from Figure 13 that the shut-in operation could only affect the pressure-temperature (P-T) changing behaviors in a short range. Therefore, we infer from the combination results from Figures 7 and 13 that the shut-in operation could hardly affect the fluid production behaviors (represented by IGWR) and formation P-T response in the long run.
All the simulation cases in this paper were conducted under the constant confining pressure and constant confining temperature conditions at lab scale. The hydrate production processes are compared qualitatively between lab scale and field scale, although the lab-scale results might be influenced by the boundary effect. The conclusions obtained in this paper would have significant importance for field hydrate production operation.

Summary and Conclusions
Both experimental and numerical simulations were conducted in the lab-scale to determine the gas-water production behaviors and formation physical responses to depressurization modes. Two experiments were used to verify the numerical model and to get the accurate value of hydrate intrinsic dissociation constant under constant confining temperature and pressure conditions. Based on the verified lab-scale numerical model, we carried out a series of numerical simulations to determine the influences of depressurization modes on cumulative gas-water production behaviors, formation temperature, and hydrate saturation. All simulation cases were conducted under the same total pressure drop (4 MPa). The influences of shut-in operation were also discussed in this paper. This study yields the following results and conclusions.
(1) The gas production behaviors match well with those of numerical simulation. However, the pressure evolutionary behaviors from experiments showed obvious hysteresis, although the trend is similar with that from numerical results. We inferred that the inlet of pressure sensor in the experiments could be blocked by secondary hydrate. The modified hydrate intrinsic dissociation constant K D is 2:43 × 10 4 (mol/m 2 ·Pa) according to the matching (2) The influences of depressurization mode on multiphase seepage field and the residual fluid (gas or water) saturation field could be characterized by the integrated gas-water ration (IGWR). The IGWR increases greatly with the increase in initial hydrate saturation. Under the same initial hydrate saturation and the same total pressure-drop magnitude, the IGWR decreases linearly with the increase in depressurizing magnitude in each step, while increase in a logarithmic style with the increase in the number of steps. Therefore, a multistep depressurization is vital for increasing the IGWR (3) Both the temperature decrease and hydrate dissociation initiate at the outlet and propagate towards the far side of the model in a nonpiston style in the I-direction. In the R-direction, the hydrate saturation and temperature evolve axisymmetrically under the combination of endothermic hydrate dissociation process, Joule-Thomson effect, and heat supply from the confining oil. Slow and multistage depressurization strategy would help to alleviate temperature drop; hence, it probably benefits for avoiding potential hydrate reformation (4) The pore pressure rebounds logarithmically if the gas production is suspended. Once the shut-in operation is ended, the pore pressure would decrease to the regular level instantaneously. Furthermore, the shut-in operation could only affect the IGWR and formation P-T response in a short range, rather than in the long-term level. However, the shut-in operation would prolong the hydrate dissociation duration greatly

Data Availability
Data are available on request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.