Numerical Simulation of Gas Hydrate Production Using the Cyclic Depressurization Method in the Ulleung Basin of the Korea East Sea

The depressurization method is known as the most productive and effective method for successful methane recovery from hydrate deposits. However, this method can cause considerable subsidence because of the increased effective stress. Maintenance of geomechanical stability is necessary for sustainable production of gas from gas hydrate deposits. In this study, the cyclic depressurization method, which uses changing the bottomhole pressure and production time during primary and secondary depressurization stage, was utilized in order to increase stability in the Ulleung Basin of the Korea East Sea. Various case studies were conducted with alternating bottomhole pressure and production time of the primary and secondary depressurization stages over 400 days. Geomechanical stability was significantly enhanced, while cumulative gas production was relatively less reduced or nearly maintained. Specially, the cumulative gas production of the 6 MPa case was more than three times higher than that of the 9 MPa case, while vertical displacement was similar between them. Therefore, it was found that the cyclic depressurization method should be applied for the sake of geomechanical stability.


Introduction
Natural gas hydrates are non-stoichiometric ice-like compounds formed by the trapping of a guest molecule (usually methane gases) in a cage of hydrogen-bonded water molecules and were discovered in 1810 by Sir Humphrey Davy [1]. Methane gas hydrate is considered to be a clean energy, because it burns cleanly and produces less CO 2 than petroleum and coal [2]. Accordingly, it can replace conventional energy sources, such as coal, petroleum and natural gas. Generally, gas hydrates are discovered in permafrost, marine sediments and deep lakes, where several necessary conditions are satisfied at the same time [3]. These conditions include low temperature (typically less than 300 K), high pressure (typically over 0.6 MPa) and the presence of methane gas and free water. Figure 1 represents that where gas hydrate has been recovered, where it is inferred to be present on the basis of seismic data and where drilling expeditions have been completed in permafrost and deep marine environments [4]. The worldwide resource of natural gas hydrate is estimated to be 3000 trillion cubic meters (TCM), which is much larger than that of conventional gas (~404 TCM) and shale gas (204~456 TCM) [5]. Marine sediments contain nearly 99 percent of global gas hydrate resources and have great potential as a new energy resources. Methane is not currently commercially produced from gas hydrates anywhere in the world; some countries have conducted field trials, for example, Messoyakha gas field in Russia, Malik site in Canada, Alaska's North Slope in USA, Nankai Trough in Japan and Northern South China Sea in China [6]. As shown in Figure 2, production technologies of the gas hydrate are related to hydrate dissociation in the P-T diagram [7]. There are three gas hydrate production methods using hydrate dissociation, namely the depressurization method, thermal stimulation method and gas-swapping or inhibitor-injection method. In the depressurization method, reservoir pressure is reduced to dissociate the gas hydrate below the phase boundary. This method is the simplest and most economical when compared to the other methods. The reservoir temperature is increased to destabilize the gas hydrate by injecting hot fluids or direct electrical heating in the thermal stimulation method. However, it is economically unfavorable, because the energy amount to increase the temperature of hydratebearing sediments (HBSs) is extremely high. In the gas-swapping or inhibitor-injection method, new hydrate phase equilibrium (dash line) is acquired by the injection of gas or inhibitor, as illustrated in Figure 2. Methane gas in hydrate is replaced with injected gas, and it is released. For example, CO2 can be used to displace methane from the hydrate while being trapped as CO2 hydrate. However, the CO2 hydrate formation can result in pore-space clogging, followed by a substantial decrease of porosity and permeability [8]. This phenomenon leads to severe loss of hydraulic conductivity. Finally, injecting inhibitors, such as methanol or brine, also dissociate hydrate. However, this method is not widely used in real cases because of non-economic and non-environmental drawbacks [9,10]. Thus, depressurization method is the best method for successful methane recovery from hydrate deposits [11,12]. As shown in Figure 2, production technologies of the gas hydrate are related to hydrate dissociation in the P-T diagram [7]. There are three gas hydrate production methods using hydrate dissociation, namely the depressurization method, thermal stimulation method and gas-swapping or inhibitor-injection method. In the depressurization method, reservoir pressure is reduced to dissociate the gas hydrate below the phase boundary. This method is the simplest and most economical when compared to the other methods. The reservoir temperature is increased to destabilize the gas hydrate by injecting hot fluids or direct electrical heating in the thermal stimulation method. However, it is economically unfavorable, because the energy amount to increase the temperature of hydrate-bearing sediments (HBSs) is extremely high. In the gas-swapping or inhibitor-injection method, new hydrate phase equilibrium (dash line) is acquired by the injection of gas or inhibitor, as illustrated in Figure 2. Methane gas in hydrate is replaced with injected gas, and it is released. For example, CO 2 can be used to displace methane from the hydrate while being trapped as CO 2 hydrate. However, the CO 2 hydrate formation can result in porespace clogging, followed by a substantial decrease of porosity and permeability [8]. This phenomenon leads to severe loss of hydraulic conductivity. Finally, injecting inhibitors, such as methanol or brine, also dissociate hydrate. However, this method is not widely used in real cases because of non-economic and non-environmental drawbacks [9,10]. Thus, depressurization method is the best method for successful methane recovery from hydrate deposits [11,12]. Messoyakha gas field in Russia, Malik site in Canada, Alaska's North Slope in USA, Nankai Trough in Japan and Northern South China Sea in China [6]. As shown in Figure 2, production technologies of the gas hydrate are related to hydrate dissociation in the P-T diagram [7]. There are three gas hydrate production methods using hydrate dissociation, namely the depressurization method, thermal stimulation method and gas-swapping or inhibitor-injection method. In the depressurization method, reservoir pressure is reduced to dissociate the gas hydrate below the phase boundary. This method is the simplest and most economical when compared to the other methods. The reservoir temperature is increased to destabilize the gas hydrate by injecting hot fluids or direct electrical heating in the thermal stimulation method. However, it is economically unfavorable, because the energy amount to increase the temperature of hydratebearing sediments (HBSs) is extremely high. In the gas-swapping or inhibitor-injection method, new hydrate phase equilibrium (dash line) is acquired by the injection of gas or inhibitor, as illustrated in Figure 2. Methane gas in hydrate is replaced with injected gas, and it is released. For example, CO2 can be used to displace methane from the hydrate while being trapped as CO2 hydrate. However, the CO2 hydrate formation can result in pore-space clogging, followed by a substantial decrease of porosity and permeability [8]. This phenomenon leads to severe loss of hydraulic conductivity. Finally, injecting inhibitors, such as methanol or brine, also dissociate hydrate. However, this method is not widely used in real cases because of non-economic and non-environmental drawbacks [9,10]. Thus, depressurization method is the best method for successful methane recovery from hydrate deposits [11,12]. However, most HBSs consist of unconsolidated porous layers, and subsidence occurs in unconsolidated sands when the reservoir pressure drops below a critical value [13,14]. Therefore, gas hydrate production that uses the depressurization method can lead to subsidence, due to the decreased strength and stiffness of HBS [15][16][17][18]. This subsidence may induce various geological disasters, such as sediment deformation, casing deformation and production platform collapse [19]. However, there have been no research studies for preventing subsidence in the case of gas hydrate production until now.
In this study, simulation studies were conducted by using the cyclic depressurization method for the sustainable gas hydrate production in the Ulleung Basin of the Korea East Sea. This method, which uses alternating depressurization and shut-in periods, was proposed for enhancing the recovery factor [20]. The simple depressurization method had a low recovery factor, because the sensible heat was not sufficiently supplied from overburden and underburden. However, the recovery factor from using the cyclic depressurization method was larger than that of the simple depressurization method. The reason is that gas hydrate was dissociated by the geothermal heat supply from overburden and underburden during the shut-in period. On the other hand, this study used the cyclic depressurization method to ensure geomechanically stable production, using high bottomhole pressure, in the secondary depressurization stage. Geomechanical stability is enhanced during the secondary depressurization stage. This study is novel in several ways. We analyzed the vertical displacement of the Ulleung Basin of the Korea East Sea during gas hydrate production, using cyclic depressurization method. Moreover, for our analysis of the vertical displacement, we conducted a reservoir simulation by using the logging data of UBGH2-6 in Ulleung Basin, both a permeability model and the relative permeability of field samples. Finally, we performed the sensitivity analysis of vertical displacement according to the cyclic bottomhole pressure and production time during primary depressurization and secondary depressurization, and it is meaningful in that it presented quantitative results of vertical displacement.

Geology of the Ulleung Basin and Hydrate Class
The Ulleung Basin is located in the Korean East Sea, between the Korean Peninsula and Japan. In order to develop gas hydrates as a potential energy source, geophysical surveys and geological studies of gas hydrates in the Ulleung Basin have been carried out since 1997. Moreover, two deep-drilling expeditions were completed in 2007 and 2010 [21]. Seismic profiles of the Ulleung Basin show that mass transport deposits and hemipelagic mud in the Ulleung Bain were generally interbedded with sandy turbidites [22]. Examinations of core and well-log data from the UBGH2 drill sites suggest that the UBGH2-6 location has relatively good gas-hydrate-reservoir-quality properties in terms of the individual and cumulative thicknesses of HBS with desirable textures (sandy) [23]. Water depth of UBGH2-6 is 2157 m, and seafloor temperature is 2.5 • C, with a geothermal gradient of 112 • C/km. The HBS is from 140 m below the seafloor (mbsf) to 153 mbsf.
UBGH2-6 site appears to be a hybrid of two hydrate classes, i.e., Class 2 and Class 3 [11,23,24]. Class 2 comprises an HBS overlying a mobile water zone, and Class 3 involves a HBS confined between two near-zero permeability strata, such as shales. Because of the very low permeability of the muds in the overburden, underburden and interlayers, the characteristics of a Class 3 deposit are dominant [25].

Reservoir Simulator
There are several hydrate reservoir simulators, including TOUGH+HYDRATE (LBNL), MH-21 (NIAIST), HydrateResSim (LBNL), STOMP (PNNL) and CMG STARS (Computer Modeling Group Ltd., Calgary, AB, Canada) [26]. Comparative previous studies have confirmed the adaptability of STARS, which can predict geomechanical system response [1,[27][28][29]. In this study, we conducted a numerical simulation of gas hydrate production, using STARS. Hydrate dissociation rate and reformation rate according to the pressure and temperature are calculated by equilibrium kinetics in STARS. Phase change from hydrate to gas is possible when using this model. Geomechanics module is a dual grid system; geomechanical and reservoir grids are independent. One-way and two-way coupling are both possible in STARS. The former does not update pore compressibility and absolute permeability after calculating displacement, stress and strain ( Figure 3). For the realistic geomechanical simulation, the two-way coupling option was used to consider updated pore compressibility and absolute permeability after geomechanical simulation.
(Computer Modeling Group Ltd., Calgary, AB, Canada) [26]. Comparative previous studies have confirmed the adaptability of STARS, which can predict geomechanical system response [1,[27][28][29]. In this study, we conducted a numerical simulation of gas hydrate production, using STARS. Hydrate dissociation rate and reformation rate according to the pressure and temperature are calculated by equilibrium kinetics in STARS. Phase change from hydrate to gas is possible when using this model. Geomechanics module is a dual grid system; geomechanical and reservoir grids are independent. One-way and two-way coupling are both possible in STARS. The former does not update pore compressibility and absolute permeability after calculating displacement, stress and strain ( Figure 3). For the realistic geomechanical simulation, the two-way coupling option was used to consider updated pore compressibility and absolute permeability after geomechanical simulation.  Figure 4 shows a schematic of the reservoir model. The cylindrical coordinated system was used to model the gas hydrate layer with overburden and underburden. Because of the generally low permeability of the system and the relative abundance of clays, well spacing was restricted to 500 m, resulting in the 250 m outer radius of the cylindrical system [23]. The overburden thickness was 140 m, and the HBS thickness was 13 m. Underburden thickness was 300 m. In order to increase numerical convergence stability, discretization along the r-direction was not uniform, instead increasing logarithmically from the center of the system. The z-direction grid size was 10 m for the overburden and underburden, and 0.1 m for the HBS. Therefore, the system was discretized into 160 × 1 × 90 = 14,400 gridblocks. The vertical production well was located at the center of the cylindrical model, and the radius of the well was 0.1 m. The geomechanical grid system was the same as the reservoir grid system.  Figure 4 shows a schematic of the reservoir model. The cylindrical coordinated system was used to model the gas hydrate layer with overburden and underburden. Because of the generally low permeability of the system and the relative abundance of clays, well spacing was restricted to 500 m, resulting in the 250 m outer radius of the cylindrical system [23]. The overburden thickness was 140 m, and the HBS thickness was 13 m. Underburden thickness was 300 m. In order to increase numerical convergence stability, discretization along the r-direction was not uniform, instead increasing logarithmically from the center of the system. The z-direction grid size was 10 m for the overburden and underburden, and 0.1 m for the HBS. Therefore, the system was discretized into 160 × 1 × 90 = 14,400 gridblocks. The vertical production well was located at the center of the cylindrical model, and the radius of the well was 0.1 m. The geomechanical grid system was the same as the reservoir grid system.  Various previous simulation studies have utilized many assumptions, such as the simple reservoir model, correlation relative permeability curve model and permeability model [12,18,23,25,[30][31][32]. Accordingly, the reservoir model of this study was constructed by using realistic field data, as shown in Figure 5 and Table 1. This model used the logging data of UBGH2-6, and the HBS of this model is composed of 14 sand layers and 13 shale layers [33]. The sand porosity is 45%, and mud porosity is 67%. The range of hydrate sat- Various previous simulation studies have utilized many assumptions, such as the simple reservoir model, correlation relative permeability curve model and permeability Appl. Sci. 2021, 11, 9748 5 of 15 model [12,18,23,25,[30][31][32]. Accordingly, the reservoir model of this study was constructed by using realistic field data, as shown in Figure 5 and Table 1. This model used the logging data of UBGH2-6, and the HBS of this model is composed of 14 sand layers and 13 shale layers [33]. The sand porosity is 45%, and mud porosity is 67%. The range of hydrate saturation is 38.8~86.2%. Moreover, we adopted experimental data of the relative permeability curve and permeability model, using a core sample from UBGH2-6, as illustrated in Figures 6 and 7 [33,34]. The experimental results of relative permeability were validated with results of X-ray CT (Computerized Tomography), and it represented good matching results [33]. Furthermore, even though the intrinsic permeability was different with each soil specimen, the permeability reduction trends with increasing hydrate saturation were similar for all samples, and the N values of Figure 7 represent the porosity [34]. Various previous simulation studies have utilized many assumptions, such as the simple reservoir model, correlation relative permeability curve model and permeability model [12,18,23,25,[30][31][32]. Accordingly, the reservoir model of this study was constructed by using realistic field data, as shown in Figure 5 and Table 1. This model used the logging data of UBGH2-6, and the HBS of this model is composed of 14 sand layers and 13 shale layers [33]. The sand porosity is 45%, and mud porosity is 67%. The range of hydrate saturation is 38.8~86.2%. Moreover, we adopted experimental data of the relative permeability curve and permeability model, using a core sample from UBGH2-6, as illustrated in Figures 6 and 7 [33,34]. The experimental results of relative permeability were validated with results of X-ray CT (Computerized Tomography), and it represented good matching results [33]. Furthermore, even though the intrinsic permeability was different with each soil specimen, the permeability reduction trends with increasing hydrate saturation were similar for all samples, and the N values of Figure 7 represent the porosity [34].     Table 1. Initial conditions and properties.

Parameter Value
Overburden thickness (m) 140 Underburden thickness (m) 300 Layer thicknesses and porosities As in Figure 5 Hydrate saturation in HBL As in Figure 5 Initial pressure at top layer (MPa) 22

Validation of the Geomechanical Model
A number of geomechanical simulation studies have been conducted for UBGH2-6 site [12,18,25,31,32] (Table 2). These studies used a different simulator, such as TOUGH+HYDRATE with FLAC2D or FLAC3D. The observation point, production time

Validation of the Geomechanical Model
A number of geomechanical simulation studies have been conducted for UBGH2-6 site [12,18,25,31,32] (Table 2). These studies used a different simulator, such as TOUGH+HYDRATE with FLAC2D or FLAC3D. The observation point, production time

Validation of the Geomechanical Model
A number of geomechanical simulation studies have been conducted for UBGH2-6 site [12,18,25,31,32] (Table 2). These studies used a different simulator, such as TOUGH+HY DRATE with FLAC2D or FLAC3D. The observation point, production time and depressurization rate were also different among studies. Our simulation results were compared with other results for general comparative purpose. One-month production simulations were conducted with various bottomhole pressure conditions, for example, 3, 6, 9 and 12 MPa. Moreover, the observation point was set at the top of HBS, and the depressurization rate was not considered.
For comparison with previous studies, the simulation results of this study are represented in Table 3 and Figure 8. The amount of vertical displacements in all cases rapidly increased and converged to a certain level ( Figure 8). As the bottomhole pressure decreased, the final amount of vertical displacement increased. The final vertical displacement was between −0.47 and −2.30 m, depending on the bottomhole pressure ( Table 3). The results of Lim [31], Kim et al. [12] and Kim [32] could not be compared with this study because of the different bottomhole-pressure conditions. From the results of Lim [31], we see that the vertical displacement was between −0.36 and −0.67 m; the results of Kim et al. [12] showed −0.4 m as a the vertical displacement. For the results of Kim [32], the range of vertical displacement was between −0.78 and −0.85 m. In this study, our simulation results were compared with the results of Moridis et al. [25] and Moon et al. [18]. The results over 14 days of production time were similar to those of Moridis et al. [25] and Moon et al. [18]. The vertical displacement value was between −0.36 and −2.06 m in this study, compared to between −0.8 and −1.6 m and between −0.63 and −1.20 m in Moridis et al. and Moon et al., respectively. In particular, the value of this study was −0.85 m, and that of Moon's study was between −0.85 and −0.89 m in the case of bottomhole pressure having 9 MPa. Therefore, we confirmed that our reservoir model and CMG reservoir simulator can appropriately predict geomechanical response.

Cyclic Depressurization Simulation Study
While the cyclic depressurization method of Konno's research consisted of a primary depressurization stage and shut-in stage [20], we used both a primary and secondary depressurization stage to produce more gas and increase pore pressure during the secondary depressurization stage. The parameters of the case study are the bottomhole pressure and production time of both the primary and secondary production stage, as shown in Table 4. The base case assumed that the bottomhole pressure and production time of the primary depressurization stage were 9 MPa and 8 days, respectively, compared to 16 MPa and 2 days for secondary depressurization stage. The shut-in of gas production was also considered solely during the secondary depressurization stage. In order to analyze an effect of the cyclic depressurization method, results of the non-cyclic depressurization case were also represented, and the production time was set as 400 days. The non-cyclic depressurization case used the values of the base case at both the primary and secondary depressurization stage.

Cyclic Depressurization Simulation Study
While the cyclic depressurization method of Konno's research consisted of a primary depressurization stage and shut-in stage [20], we used both a primary and secondary depressurization stage to produce more gas and increase pore pressure during the secondary depressurization stage. The parameters of the case study are the bottomhole pressure and production time of both the primary and secondary production stage, as shown in Table 4. The base case assumed that the bottomhole pressure and production time of the primary depressurization stage were 9 MPa and 8 days, respectively, compared to 16 MPa and 2 days for secondary depressurization stage. The shut-in of gas production was also considered solely during the secondary depressurization stage. In order to analyze an effect of the cyclic depressurization method, results of the non-cyclic depressurization case were also represented, and the production time was set as 400 days. The non-cyclic depressurization case used the values of the base case at both the primary and secondary depressurization stage.

No. Stage Property Value
Case study 3 (Figures 13 and 14

Bottomhole Pressure during Primary Depressurization Stage
The bottomhole pressure during the primary depressurization was varied (6, 9 and 12 MPa). The case of the lower bottomhole pressure caused a high gas production rate and cumulative gas production ( Figure 9). The gas-production rate of all cases steadily increased until 60 days, and it maintained a certain level, with the exclusion of the 6 MPa case (Figure 9a). The gas-production rate of the 6 MPa case consistently decreased after 60 days, and that of the 12 MPa case was very low compared to other cases. The cumulative gas production was between 1.64 × 10 4 m 3 (in the case of 12 MPa) and 7.47 × 10 5 m 3 (in the case of 6 MPa) (Figure 9b). There was little difference in the cumulative gas production between the non-cyclic case and the 9 MPa case; the difference was only 12.1%. Specially, the cumulative gas production of the 6 MPa case was three times higher than that of the 9 MPa case. During the secondary pressure stage, a small amount of gas production cyclically decreased in all cases because of the relatively high bottomhole pressure condition.

Bottomhole Pressure during Primary Depressurization Stage
The bottomhole pressure during the primary depressurization was varied (6, 9 and 12 MPa). The case of the lower bottomhole pressure caused a high gas production rate and cumulative gas production ( Figure 9). The gas-production rate of all cases steadily increased until 60 days, and it maintained a certain level, with the exclusion of the 6 MPa case (Figure 9a). The gas-production rate of the 6 MPa case consistently decreased after 60 days, and that of the 12 MPa case was very low compared to other cases. The cumulative gas production was between 1.64 × 10 4 m 3 (in the case of 12 MPa) and 7.47 × 10 5 m 3 (in the case of 6 MPa) (Figure 9b). There was little difference in the cumulative gas production between the non-cyclic case and the 9 MPa case; the difference was only 12.1%. Specially, the cumulative gas production of the 6 MPa case was three times higher than that of the 9 MPa case. During the secondary pressure stage, a small amount of gas production cyclically decreased in all cases because of the relatively high bottomhole pressure condition.
(a) (b) Figure 9. Results of gas production by use of different bottomhole pressure during primary depressurization stage: (a) gas production rate and (b) cumulative gas production. Figure 10 represents the vertical subsidence at the top of the HBS. The amount of vertical displacement increased according to the production time. In the case of low bottomhole pressure, the amount of vertical displacement was high due to the low pore pressure; the range of vertical displacement was from −1.09 m (in the case of 12 MPa) to −2.39 m (in the case of 6 MPa). From the aforementioned cumulative gas-production results, we confirmed that high cumulative gas production resulted in a high amount of vertical subsidence. In all cases, the amount of vertical displacement increased during the primary depressurization stage, while it decreased during the secondary depressurization stage. The reason is that relatively low gas production during the secondary depressurization stage caused the increment of pore pressure as compared to the primary depressurization stage. In addition, the amount of vertical subsidence in the case of 9 MPa was low com- Figure 9. Results of gas production by use of different bottomhole pressure during primary depressurization stage: (a) gas production rate and (b) cumulative gas production. Figure 10 represents the vertical subsidence at the top of the HBS. The amount of vertical displacement increased according to the production time. In the case of low bottomhole pressure, the amount of vertical displacement was high due to the low pore pressure; the range of vertical displacement was from −1.09 m (in the case of 12 MPa) to −2.39 m (in the case of 6 MPa). From the aforementioned cumulative gas-production results, we confirmed that high cumulative gas production resulted in a high amount of vertical subsidence. In all cases, the amount of vertical displacement increased during the primary depressurization stage, while it decreased during the secondary depressurization stage. The reason is that relatively low gas production during the secondary depressurization stage caused the increment of pore pressure as compared to the primary depressurization stage. In addition, the amount of vertical subsidence in the case of 9 MPa was low compared to that of the non-cyclic case; the difference was only 16.6%. This value was higher than the difference of the cumulative gas production. Moreover, the result of the 6 MPa case was similar with that of the non-cyclic case. Accordingly, geomechanical stability increased significantly by using the cyclic depressurization method, and this is a key parameter for predicting geomechanical stability.

Results of Production Time Case during Primary Depressurization Stage
The production time during the primary depressurization stage ranged from 2 to 8 days. As represented in Figure 11a, more gas was produced in all cyclic depressurization cases than the non-cyclic case during the primary depressurization stage, and the gas production rate of all cases was kept at a certain level. At the initial production time, the gas production rate during the secondary depressurization stage increased in all cases. The maximum gas production rate was about 193 m 3 /day in the case of 2 days. However, the gas production rate reduced after 100 days during the secondary depressurization stage. As production days of primary depressurization increased, cumulative gas production also increased from 1.43 × 10 5 m 3 (in the case of 2 days) to 2.00 × 10 5 m 3 (in the case of 8 days) (Figure 11b). The cumulative gas production of the non-cyclic case was higher than that of the cyclic cases. In the case of 2 days, cumulative gas production was less than half that of the non-cyclic case.
(a) (b) Figure 11. Results of gas production by use of different production time during primary depressurization stage: (a) gas production rate and (b) cumulative gas production.

Results of Production Time Case during Primary Depressurization Stage
The production time during the primary depressurization stage ranged from 2 to 8 days. As represented in Figure 11a, more gas was produced in all cyclic depressurization cases than the non-cyclic case during the primary depressurization stage, and the gas production rate of all cases was kept at a certain level. At the initial production time, the gas production rate during the secondary depressurization stage increased in all cases. The maximum gas production rate was about 193 m 3 /day in the case of 2 days. However, the gas production rate reduced after 100 days during the secondary depressurization stage. As production days of primary depressurization increased, cumulative gas production also increased from 1.43 × 10 5 m 3 (in the case of 2 days) to 2.00 × 10 5 m 3 (in the case of 8 days) (Figure 11b). The cumulative gas production of the non-cyclic case was higher than that of the cyclic cases. In the case of 2 days, cumulative gas production was less than half that of the non-cyclic case.
During the secondary depressurization stage, the vertical subsidence of the cyclic cases decreased more than the non-cyclic case, due to the increased pore pressure, for example, from −2.29 m (in the case of non-cyclic) to between −1.56 m (in the case of 2 days) and −1.91 m (in the case of 8 days) after 400 days ( Figure 12). Moreover, our simulation results showed that short-production-time cases (i.e., 2-day case) provided good stability, showing the low vertical subsidence due to the low gas production. maximum gas production rate was about 193 m /day in the case of 2 days. However, the gas production rate reduced after 100 days during the secondary depressurization stage. As production days of primary depressurization increased, cumulative gas production also increased from 1.43 × 10 5 m 3 (in the case of 2 days) to 2.00 × 10 5 m 3 (in the case of 8 days) (Figure 11b). The cumulative gas production of the non-cyclic case was higher than that of the cyclic cases. In the case of 2 days, cumulative gas production was less than half that of the non-cyclic case.
(a) (b) Figure 11. Results of gas production by use of different production time during primary depressurization stage: (a) gas production rate and (b) cumulative gas production.
During the secondary depressurization stage, the vertical subsidence of the cyclic cases decreased more than the non-cyclic case, due to the increased pore pressure, for example, from −2.29 m (in the case of non-cyclic) to between −1.56 m (in the case of 2 days) Figure 11. Results of gas production by use of different production time during primary depressurization stage: (a) gas production rate and (b) cumulative gas production. and −1.91 m (in the case of 8 days) after 400 days ( Figure 12). Moreover, our simulation results showed that short-production-time cases (i.e., 2-day case) provided good stability, showing the low vertical subsidence due to the low gas production.

Results of Bottomhole Pressure Case during Secondary Depressurization Stage
The bottomhole pressures during secondary depressurization stage varied from 12 to 20 MPa, except for in the shut-in case. The gas-production rates of all cases were similar and maintained a certain value (Figure 13a). In the case of the shut-in and 20 MPa, gasproduction rates during the secondary depressurization stage were nearly zero. On the other hand, the 12 MPa case showed that the high gas-production rate and maximum gasproduction rate were 385 m 3 /day after 110 days. Then, the gas production rate reduced, and it was still 165 m 3 /day after 400 days. There were only small differences in the cumulative gas production (Figure 13b). For example, each cumulative gas production was between 2.00 × 10 5 m 3 and 2.15 × 10 5 m 3 , while that was 2.28 × 10 5 m 3 for the non-cycle case; herein, the cumulative gas production of the 12 MPa case was 2.15 × 10 5 m 3 compared to 2.00 × 10 5 m 3 for the 16 MPa case (base case). It increased by 7.3% due to the high gasproduction rate during the secondary depressurization stage. Specially, an interesting fact is that the cumulative gas production of both the shut-in case and 20 MPa case was higher than that of the 16 MPa case, although no gas was produced in the case of the shut-in, and little gas was produced in the 20 MPa case during the secondary stage. The reason may be that more geothermal heat was supplied from overburden and underburden during the shut-in period.

Results of Bottomhole Pressure Case during Secondary Depressurization Stage
The bottomhole pressures during secondary depressurization stage varied from 12 to 20 MPa, except for in the shut-in case. The gas-production rates of all cases were similar and maintained a certain value (Figure 13a). In the case of the shut-in and 20 MPa, gasproduction rates during the secondary depressurization stage were nearly zero. On the other hand, the 12 MPa case showed that the high gas-production rate and maximum gasproduction rate were 385 m 3 /day after 110 days. Then, the gas production rate reduced, and it was still 165 m 3 /day after 400 days. There were only small differences in the cumulative gas production (Figure 13b). For example, each cumulative gas production was between 2.00 × 10 5 m 3 and 2.15 × 10 5 m 3 , while that was 2.28 × 10 5 m 3 for the non-cycle case; herein, the cumulative gas production of the 12 MPa case was 2.15 × 10 5 m 3 compared to 2.00 × 10 5 m 3 for the 16 MPa case (base case). It increased by 7.3% due to the high gas-production rate during the secondary depressurization stage. Specially, an interesting fact is that the cumulative gas production of both the shut-in case and 20 MPa case was higher than that of the 16 MPa case, although no gas was produced in the case of the shut-in, and little gas was produced in the 20 MPa case during the secondary stage. The reason may be that more geothermal heat was supplied from overburden and underburden during the shut-in period. Appl. Sci. 2021, 11, x FOR PEER REVIEW 12 of 15 (a) (b) Figure 13. Results of gas production by use of different bottomhole pressure during secondary depressurization stage: (a) gas production rate and (b) cumulative gas rate.
The vertical displacement of bottomhole pressure during secondary depressurization stage was between −1.56 m (shut-in case) and −1.91 m (in the case of 12 MPa) after 400 days, and the geomechanical stability in all cases increased ( Figure 14). This parameter had little effect on vertical displacement during the primary depressurization stage according to the bottomhole pressure, while the vertical displacement strongly restored during secondary depressurization stage in the case of high bottomhole pressure (16 and 20 MPa) or shut-in. The reason is that little gas was produced in the case of the high bottomhole pressure. In terms of gas production, the low-bottomhole-pressure case was more productive, although the geomechanical stability was not good.

Results of Production Time Case during Secondary Depressurization Stage
Simulations were conducted with production days during the secondary depressurization stage changed by 1, 2 and 4 days. As shown in Figure 15a, the production rates during the primary depressurization stage were similar for all cases. In the case of 1 day, maximum gas production rate was about 126 m 3 /day and there was little gas production after 100 days during secondary depressurization stage. The gas-production rate of the 4day case during the secondary depressurization stage was high at the initial production The vertical displacement of bottomhole pressure during secondary depressurization stage was between −1.56 m (shut-in case) and −1.91 m (in the case of 12 MPa) after 400 days, and the geomechanical stability in all cases increased ( Figure 14). This parameter had little effect on vertical displacement during the primary depressurization stage according to the bottomhole pressure, while the vertical displacement strongly restored during secondary depressurization stage in the case of high bottomhole pressure (16 and 20 MPa) or shut-in. The reason is that little gas was produced in the case of the high bottomhole pressure. In terms of gas production, the low-bottomhole-pressure case was more productive, although the geomechanical stability was not good. The vertical displacement of bottomhole pressure during secondary depressurization stage was between −1.56 m (shut-in case) and −1.91 m (in the case of 12 MPa) after 400 days, and the geomechanical stability in all cases increased ( Figure 14). This parameter had little effect on vertical displacement during the primary depressurization stage according to the bottomhole pressure, while the vertical displacement strongly restored during secondary depressurization stage in the case of high bottomhole pressure (16 and 20 MPa) or shut-in. The reason is that little gas was produced in the case of the high bottomhole pressure. In terms of gas production, the low-bottomhole-pressure case was more productive, although the geomechanical stability was not good.

Results of Production Time Case during Secondary Depressurization Stage
Simulations were conducted with production days during the secondary depressurization stage changed by 1, 2 and 4 days. As shown in Figure 15a, the production rates during the primary depressurization stage were similar for all cases. In the case of 1 day, maximum gas production rate was about 126 m 3 /day and there was little gas production after 100 days during secondary depressurization stage. The gas-production rate of the 4day case during the secondary depressurization stage was high at the initial production

Results of Production Time Case during Secondary Depressurization Stage
Simulations were conducted with production days during the secondary depressurization stage changed by 1, 2 and 4 days. As shown in Figure 15a, the production rates during the primary depressurization stage were similar for all cases. In the case of 1 day, maximum gas production rate was about 126 m 3 /day and there was little gas production after 100 days during secondary depressurization stage. The gas-production rate of the 4-day case during the secondary depressurization stage was high at the initial production time, and the maximum gas production rate was 214 m 3 /day; however, it reduced after 96 days. Each cumulative gas production was between 1.81 × 10 5 m 3 and 2.27 × 10 5 m 3 ( Figure 15b). Specially, cumulative gas production of the 1-day case was very similar to the production of the non-cyclic case, and their difference was only 0.5%. time, and the maximum gas production rate was 214 m 3 /day; however, it reduced after 96 days. Each cumulative gas production was between 1.81 × 10 5 m 3 and 2.27 × 10 5 m 3 ( Figure  15b). Specially, cumulative gas production of the 1-day case was very similar to the production of the non-cyclic case, and their difference was only 0.5%.
(a) (b) Figure 15. Results of gas production by use of different production time during secondary depressurization stage. (a) Gas production rate; (b) Cumulative gas rate.
The vertical subsidence results are represented in Figure 16. The final vertical displacement values during primary depressurization stage were between −1.87 m (in the case of 4 days) and −1.92 m (in the case of 1 day), and the difference was very small. However, the final vertical displacement during the secondary depressurization stage was between −1.19 m (in the case of 4 days) and −1.65 m (in the case of 1 day), and the difference was large. The reason is that pore pressure was restored consistently because the low gas production time was longer. Therefore, additional production time in the secondary stage did not increase stability, because the final vertical displacement values were similar during primary depressurization stage.  The vertical subsidence results are represented in Figure 16. The final vertical displacement values during primary depressurization stage were between −1.87 m (in the case of 4 days) and −1.92 m (in the case of 1 day), and the difference was very small. However, the final vertical displacement during the secondary depressurization stage was between −1.19 m (in the case of 4 days) and −1.65 m (in the case of 1 day), and the difference was large. The reason is that pore pressure was restored consistently because the low gas production time was longer. Therefore, additional production time in the secondary stage did not increase stability, because the final vertical displacement values were similar during primary depressurization stage.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 15 time, and the maximum gas production rate was 214 m 3 /day; however, it reduced after 96 days. Each cumulative gas production was between 1.81 × 10 5 m 3 and 2.27 × 10 5 m 3 ( Figure  15b). Specially, cumulative gas production of the 1-day case was very similar to the production of the non-cyclic case, and their difference was only 0.5%.
(a) (b) Figure 15. Results of gas production by use of different production time during secondary depressurization stage. (a) Gas production rate; (b) Cumulative gas rate.
The vertical subsidence results are represented in Figure 16. The final vertical displacement values during primary depressurization stage were between −1.87 m (in the case of 4 days) and −1.92 m (in the case of 1 day), and the difference was very small. However, the final vertical displacement during the secondary depressurization stage was between −1.19 m (in the case of 4 days) and −1.65 m (in the case of 1 day), and the difference was large. The reason is that pore pressure was restored consistently because the low gas production time was longer. Therefore, additional production time in the secondary stage did not increase stability, because the final vertical displacement values were similar during primary depressurization stage.

Conclusions
In this study, a field-scale numerical simulation study was conducted by using the cyclic depressurization method for the sustainable gas hydrate production. Geological model of UBGH2-6 site was constructed based on the input data from various papers and experimental data. STARS of CMG was utilized as a reservoir simulator that can analyze fluid flow, heat transfer and hydrate dissociation behavior. The geomechanical results were in good agreement with previous studies. Case studies were conducted according to bottomhole pressure and production time during both the primary and secondary depressurization stages. Geomechanical stability was enhanced during the secondary depressurization stage, and amount of vertical displacement of each case was different respectively. Specially, the 6 MPa case during primary depressurization case showed that vertical displacement was similar with the non-cyclic case, while the cumulative gas production of the former was more than three times higher than that of the latter. In addition, in the case of 1 day during the secondary depressurization stage, the cumulative gas production was almost the same as the no-cyclic case, but the geomechanical stability was more enhanced than the non-cyclic case. Accordingly, the geomechanical stability was acquired by using the cyclic depressurization method with little loss of gas production. Our next research will be conducted with a non-cylindrical field-scale reservoir model, considering multiple wells, using the cyclic depressurization method.