Numerical Evaluation of a Novel Development Mode for Challenging Oceanic Gas Hydrates Considering Methane Leakage

: The exploitation of challenging oceanic gas hydrate reservoirs with low permeability and permeable boundary layers faces the challenges of methane leakage and low production. Considering this aspect, a novel ﬁve-spot injection–production system combined with hydraulic fracturing was proposed. In particular, the potential of this development mode, including hydrate dissociation, gas production, and gas capture, was evaluated in comparison with a three-spot injection–production system. The results showed that increasing the fracture conductivity cannot prevent CH 4 leakage in the three-spot, and the leakage accounted for 5.6% of the total gas production, even at the maximum fracture conductivity of 40 D · cm. Additionally, the leakage amount increased as the well spacing increased, and the leakage accounted for 36.7% of the total gas production when the well spacing was 140 m. However, the proposed development mode completely addressed CH 4 leakage and signiﬁcantly increased gas production. The average gas production rate reached 142 m 3 /d per unit length of the horizontal section, which was expected to reach the commercial threshold. The variance analysis indicated that optimal plans for the challenging hydrates in the Shenhu area were well spacing of 100–120 m and fracture conductivity greater than 20 D · cm.


Introduction
Natural gas hydrate (NGH) is a cage-structured crystalline solid formed by natural gas (primarily methane) and water under high pressure and low temperature [1]. The global natural gas reserves in NGH are approximately 2.5-120 × 10 15 stm 3 , nearly twice the total carbon content of fossil fuels [2][3][4]; thus, NGH is generally considered to play a crucial role in future energy supply and climate management [5,6].
NGH is contended to be a challenging natural gas reservoir because its development faces complex phase transition, heat transfer, multiphase flow, formation deformation, and many other challenges, compared with conventional natural gas development [6][7][8][9][10]. According to the reservoir properties and occurrence environment, the challenging NGH reservoirs are characterized by low initial reservoir temperature (CH-k), lack of impermeable boundary layers (CH-b), high hydrate stability (CH-s), and low reservoir permeability (Ch-k) [11]. A vast majority of global oceanic NGH reservoirs have both low reservoir permeability and permeable boundary layers (CH-kb) [12,13], which is the main type of oceanic NGH and the primary target for commercial development.
The premise of efficient gas recovery from NGH reservoirs is to obtain a considerable hydrate decomposition rate. However, the depressurization method, which is generally considered to be the most promising, exhibits a low hydrate decomposition rate and poor gas production in CH-kb [11,14,15]. This is caused by two mechanisms: first, the low reservoir permeability and permeable boundaries inhibit pressure drop transmission and fluid flow [16][17][18]; second, insufficient reservoir sensible heat further reduced the hydrate decomposition capacity [19,20]. Therefore, exploring the efficient development mode of CH-kb has recently become an urgent topic.
Depressurization combined with thermal fluid injection could provide extra heat for hydrate dissociation, which is a promising improved scheme for NGH development [21,22]. In particular, the injection-production system with multiple wells, such as two-spot [23] and three-spot (TS) well patterns [24], significantly increases drainage area and enhances gas recovery, and has been envisaged by Japan for the commercial development of NGH [25]. However, its direct applicability in CH-kb is problematic because the gas released by hydrate thermal decomposition is difficult to pass through the inter-well hydrate undissociated zones where hydrates fill in pores and thereby form barriers to gas production. At present, using hydraulic fracturing (HF), a key stimulation technology in the commercial development of shale and tight gas, to create high conductivity artificial fractures between the injection well (well i ) and production well (well p ) is a promising solution. Several numerical simulations confirmed that fractures significantly enhanced hydrate decomposition and gas production capacity, wherein hydraulic fracturing combined with a three-spot well pattern (HF + TS) was used [26][27][28]. However, existing studies have not yet focused on the risk of methane leakage in hydraulic fractured CH-kb. When the hydrate decomposition front breaks through to the permeable boundaries, the gas released by hydrate thermal decomposition escapes into boundary layers under the drive of high injection pressure [29]. This is unfavorable for gas recovery and potentially pollutes the oceanic and atmospheric environment. Recently, methane released by oceanic hydrate dissociation is considered to play an important role in the global carbon cycle, as it is leaking through permeable seabed or natural fractures [30][31][32][33][34]. Furthermore, the massive release of methane induced by NGH development in the future may affect the atmospheric CH 4 budget and exacerbate the greenhouse effect [35][36][37]; therefore, there is an urgent need to evaluate the risk of methane leakage in hydraulic fractured CH-kb.
In our previous studies, we successfully addressed methane leakage in the nonfractured CH-kb using a novel five-spot well pattern, wherein two capture wells (well c ) were arranged at the boundary layers to recover the escaped gas [29]. However, because the low reservoir permeability weakens the inter-well interaction, the spacing between injection and production wells is limited to 50 m. Small space well pattern significantly increases the drilling cost, reduces economic benefits, and is not conducive to the commercial development of CH-kb. Given this, a novelty development mode combining hydraulic fracturing and the five-spot injection-production system (HF + FS) is proposed here to address the aforementioned challenges. In this study, the production potential of the proposed development mode (HF + FS) for CH-kb at the SH2 site in the Shenhu area, South China Sea, was evaluated by comparing it with HF + TS. Specifically, the impacts of fracture conductivity (F c ) and well spacing (W s ) on hydrate decomposition, gas production, methane escape, and gas recovery under the two development modes were thoroughly discussed, and ideal development plans were obtained through variance analysis (ANOVA). These findings provide valuable references for the safe and efficient development of CH-kb.

Development Plan Design
The implementation of HF + FS includes the following steps: deploying the five-spot well pattern, performing hydraulic fracturing between the well i and well p , and production with depressurization combined with thermal fluid injection.
As shown in Figure 1, all wells are horizontal and extend along the Y direction. The well i and well c are arranged by drilling a three-branch horizontal well, wherein the middle section is used for injection, and the upper and lower sections are used for gas capture.
Additionally, well i and well p are connected through a horizontal artificial fracture created by hydraulic fracturing, and its width is 10 mm [38,39]. Due to the lack of engineering cases, the injection pressure and temperature in existing simulations are [16][17][18][19][20][21][22] MPa and 30-90 • C, respectively [21][22][23][24][26][27][28][29]. For comparison, we adopted the same injection-production scheme as that in our previous study [29], that is, the injection temperature and pressure were 60 • C and 16 MPa, and the bottom hole pressures of well p , upper well c , and lower well c were 4, 10, and 11 MPa, respectively. Additionally, geomechanical issues are not considered in the model, because current research mainly focuses on CH 4 migration and hydrate decomposition behavior under the two development modes.

Generalization of the Fracture
In this study, the fracture is generalized as porous media with high conductivity, and the multi-flow in the fracture is described as a continuous flow that follows Darcy's law [40]. Therefore, the two-phase seepage in the matrix and fracture that considers the gravity and capillary forces can be given by Ref. [41]: where F l and F g i are the fluxes of the liquid and gas phases, respectively, wherein l, g, and i represent the liquid phase, gas phase, and gas component, respectively; k is the absolute permeability of the matrix, m 2 ; k l r and k g r are the relative permeabilities; µ l and µ g are the fluid viscosities, Pa·s; P l and P g are the fluid pressures, Pa; ρ l and ρ g are the fluid densities, kg/m 3 ; g is the gravitational acceleration, m/s 2 ; C g i is the mass friction; and ω is the Klinkenberg factor, Pa.
The k l r and k g r are expressed as [41]: where S l and S g are the water and gas saturation, S l ir and S g ir are the irreducible saturation of liquid and gas, n l and n g are the stone indexes of liquid and gas.
The value of P l can be obtained by Ref. [42]: where P c is the capillary pressure, Pa, and is expressed as [43]: where P t is threshold pressure, Pa; and ϑ is the pore structure index.

Numerical Code
To simulate the production behavior of CH-kb during depressurization and thermal fluid injection, the hydrate reaction based on thermodynamic conditions, mass conservation, energy conservation, and the flux of possible components and phases should be considered. The details of these governing equations have been given in previous studies [23,41]. Several simulators have been developed to solve these governing equations, such as Tough + Hydrate (T + H), and MH21-Hydrate [38]. T + H is developed by Lawrence Berkeley Laboratory and is popular in simulating field-scale hydrate reservoir development [41]. Three forms of heat exchange, heat conduction, heat convection, and heat radiation, are considered in T + H. Heat sources include hydrate reaction heat, phase change latent heat, the sensible heat of the reservoir itself, etc. The fluid flow in T + H follows Darcy's law, and the convection and molecular diffusion of dissolved gases and inhibitors in water are considered. The code includes four components, including water, CH 4 , CH 4 hydrate, and inhibitor, and four phases, including gas, liquid, ice, and hydrate. Hydrate decomposition is based on the equilibrium model or kinetic model [44]. Therefore, T + H is powerful in modeling non-isothermal hydrate decomposition, phase behavior, heat exchange, and fluid flow in hydrate reservoirs. Furthermore, the reliability of this code has been verified both at the laboratory and field scales [45,46]; thus, T + H is employed in this study.

Initial Conditions and Model Parameters
Post-drilling data at the SH2 site indicates that hydrates occur in argillaceous silt sediment below 185 mbsf (meter below seafloor), with thickness, porosity, and hydrate saturation of 0-40 m, 0.33-0.48, and 26-48%, respectively (40 m, 0.38, and 44%, respectively, in this study) [47]. Core testing revealed that the absolute permeabilities of the overburden, hydrate-bearing layer (HBL), and underburden were 10 mD [48]. Assuming that the reservoir properties are homogeneous along the Y direction, horizontal wells with long production sections are used for gas recovery and fluid injection, thus the thickness of the model in the Y direction could be set to 1 m, as shown in Figure 1.
The reservoir temperature and geothermal gradients were 14.87 • C and 0.047 • C/m, respectively [49]; therefore, the domain temperature can be obtained. The pore water is assumed to be interconnected due to the permeable boundary layers, and the pressure distribution was derived from the reservoir pressure of 14.97 MPa and pressure gradient of 0.01 MPa/m [49]. The temperature and pressure at the upper and lower boundaries of the model are considered constant because HBL is covered by seawater and underlain by enough thick sediment. Based on this, the initial equilibrium state of the reservoir is obtained ( Figure 2). The values of other model parameters are listed in Table 1.

Simulation Scenarios
The fracture conductivity depends on the stress condition, proppant performance, etc., mostly in dozens of D cm [34][35][36][37]. To obtain an ideal development plan, four levels of fracture conductivity are set at 5, 10, 20, and 40 D·cm, and four levels are set at 80, 100, 120, and 140 m for well spacing, as presented in Table 2.

Comparison of Different Development Modes
Based on the simulation results of scenarios 01, 04, and 18#, we compared the production potential of different development modes (TS, HF + TS, HF + FS) for CH-kb. In TS, the high-temperature (higher than 30 • C) zone initially spreads circularly around well i (Figure 3a,b). At 1800 d, the thermal decomposition front advanced to the boundaries (Figure 3c), and hot water flowed into the boundary layers, resulting in non-uniform heat diffusion ( Figure 3c). Compared to that in TS, the temperature evolution in HF + TS differed considerably. First, hot water mainly flowed into the fracture, resulting in a faster horizontal heat diffusion rate. Second, the hot water injection was strengthened by the fracture, leading to a larger heating area. At 360 d, the horizontal and vertical heating radii are 68 and 9 m, respectively (Figure 3d), which are higher than that at 300 d (7 m) in TS. At 600 d, the high temperature around the well p indicates that hot water has flowed into it (Figure 3e). This is a phenomenon known as water-flooding in the petroleum industry that dramatically increases the water yield, and the time of occurrence is T f . At 1800 d, the vertical heating radius reached 44 m (Figure 3f). In HF + FS, vertical heat diffused faster at 600 d ( Figure 3h). This is because the low pressure of well c induces more hot water to migrate vertically. Moreover, since the well c prevents hot water from passing through the permeable boundaries, the vertical heating radius is still 20 m at 1800 d ( Figure 2i). Therefore, the heating efficiency is significantly enhanced by the fracture, and hot water infiltration into the boundary layer is avoided by well c .  Figure 4 shows the reservoir pressure evolution in different development modes. In TS, two low-pressure (≤13 MPa) zones and a high-pressure (≥15 MPa) zone are formed around well p and well i (Figure 4a-c), respectively. As the pressure transmission and fluid injection are inhibited by the low reservoir permeability, both the low-and highpressure zones exhibit a slow diffusion rate. In HF + TS, the high-pressure zone diffuses faster along the fracture (Figure 4d). Moreover, the radius of the low-pressure zone is 48 m at 1800 d (Figure 4f), which is 2.09 times that in TS. At 1800 d, the radius of the high-pressure zone shrinks to 15 m due to the pressure relief (Figure 4f). In HF + FS, the low-pressure transmission of well p is the same as that in HF + TS, and two low-pressure zones are also formed around well c , which intersect with the low-pressure zones formed by well p (Figure 4i). This indicates that the gas is well recovered, thereby relieving the reservoir pressure.  Figure 5 shows the hydrate saturation evolution in different development modes. In TS, the depressurization decomposition zone and thermal decomposition zone expand circularly with a slow expansion rate (Figure 5a,b). At 1800 d, the decomposition radii of the thermal stimulation and depressurization are 13 and 21 m (Figure 5c), respectively, indicating that thermal stimulation has a better hydrate decomposition efficiency. Moreover, the gas leaking into the overburden recombines with free water and re-form hydrates (Figure 5c). In HF + TS, hydrates dissociate rapidly along the fracture, while more hydrates are generated in the overburden (Figure 5d-f). In HF + FS, the hydrate saturation evolution is similar to that in HF + TS, but no hydrate in the overburden (Figure 5g-i). Therefore, the arrangement of well c effectively avoids the secondary formation of hydrate in the overburden.  Figure 6 displays the gas saturation evolution in different development modes. In TS, gas saturation evolution is similar to that of hydrate saturation. The gas released by thermal stimulation accumulated at the decomposition front, and can not be produced (Figure 6a-c). Thus, there is no synergy effect between well i and well p , and thermally decomposed gas migrates into the boundary layer (Figure 6c). In HF + TS, the gas released by thermal stimulation migrates to the well p through the fracture (Figure 6d-f). However, more gas leaked into the boundary layers (Figure 6d-f), and the farthest vertical escape position was 46 m away from welli (Figure 6f), which was higher than that in TS (35 m). In HF + FS, the gas saturation evolution is similar to that in HF + TS, but there is no methane in the boundary layers (Figure 6g-i). Therefore, methane leakage is addressed by well c .  Figure 7 shows the variation of the CH 4 released rate (Q r ) and hydrate dissociation ratio (D r ) in different development modes, wherein D r is the ratio between the cumulative decomposed hydrate mass (M d , kg) and the original hydrate mass (M o , kg) in the reservoir, and is expressed as: In TS, Q r is minimum and stable at approximately 50 m 3 /d. In HF + TS, the variation of Q r exhibited three stages. In the first stage (0-50 d), the stable state of the hydrates was broken instantaneously; thus, the initial Q r value was the highest. In the second stage (50-450 d), Q r stabilized around 200 m 3 /d. In the third stage (450-1800 d), Q r gradually decreases. This is caused by two factors: first, hot water bursts into the well p along the fracture, which reduces the heating efficiency. Second, the hydrate decomposition front gradually moves away from the fracture. In HF + FS, the evolution of Q r was similar to that in HF + TS, but it was higher in the first two stages and declined faster in the third stage. This is because the presence of the well c further strengthened the initial injection, while the subsequent water-flooded further reduced the heating efficiency. The final D r in TS, HF + TS, and HF + FS is 31%, 93%, and 95%, respectively. This indicates that the final D r is improved by the fracture, whereas it is barely affected by the existence of well c . Figure 8 displays the variation of the gas production rate (Q p ) and cumulative gas production (V p ) in different development modes. In TS, Q p is minimum and totally contributed by depressurization. In HF + TS, affected by Q r , Q p also presents three stages. In the first stage (0-85 d), Q p decreases rapidly. Subsequently (85-374 d), gas migrates to well p ; thus, Q p increases and reaches its peak of 143 m 3 /d, which is 5.72 times that in TS. In the third stage (374-1800 d), the decrease of Q r leads to a decrease in Q p . For HF + FS, the variation of Q p is similar to that in HF + TS. Attributed to the gas recovery by well c , the maximum Q p reaches 203 m 3 /d, which is 1.43 times that in HF + TS. The final V p in HF + FS is 20.75 × 10 4 m 3 , which is 5.47 and 1.18 times of TS and HF + TS, respectively. Therefore, well c and the fracture synergistically enhance the production.

Methane Leakage, Gas Capture, and Recovery Efficiency
The volume of methane leakage (V l , m 3 ) consists of two parts: the volume of free gas in the boundary layers (V f l , m 3 ) and the volume of gas trapped in re-formed hydrates (V h l , m 3 ), and is expressed as: To characterize the methane capture ability of well c , the cumulative gas production of well c was defined as the volume of captured gas (V c ). Recovery efficiency (R e ) is defined as the ratio between cumulative gas production (V p , m 3 ) and the volume of gas released by hydrate decomposition (V r , m 3 ), and is expressed as: Figure 9 depicts V l in TS and HF + TS, V c in HF + FS, and R e in different development modes at 1800 d. The V l in TS and HF + TS are 15,080 and 18,389 m 3 , respectively, indicating that CH 4 leakage is aggravated by the fracture. The V c in HF + FS is 40,750 m 3 and there is no CH 4 leakage, which is consistent with the gas distribution shown in Figure 7i. Moreover, V c is higher than V l , indicating that well c also enhances V p . The R e in HF + FS is 0.93, an increase of 0.14 and 0.42 relative to TS and HF + TS, respectively. Therefore, HF + FS is a better development mode for CH-kb. It is important to note that 7% of the gas was trapped in the reservoir at 1800 d (Figure 6i). To avoid methane leakage, depressurized production can be maintained until free gas is completely produced.

The Effect of Fracture Conductivity
Scenarios 02-05, 10, 14, 18, and 22# simulated the production behavior of the two development modes (HF + TS, HF + FS) with different fracture conductivity at the well spacing of 100 m. By comparing these simulation results, the effect of fracture conductivity on production behavior is analyzed. Figure 10 depicts the variation of Q r in HF + TS and HF + FS at different F c . In HF + TS, Q r varies gently around 90 m 3 /d when F c is 5 D·cm, which is similar to that in TS, indicating that an F c value of 5 D·cm is insufficient for gas migration. When F c is 10 D·cm, T f is 1000 d. For F c of 20 and 40 D·cm, T f further increases to 400 and 70 d, respectively. Additionally, Q r at 1800 d is lower when F c ≥ 20 D·cm than that when F c ≤10 D·cm. This is attributed to the inadequacy of remaining hydrates for the gas supply when F c ≥ 20 D·cm. In HF + FS, Q r is higher at the early stage of production and lower at the later stage of production. The final D r values in HF + TS are 0.51, 0.82, 0.93, and 0.96 for F c of 5, 10, 20, and 40 D·cm, respectively (Figure 11), indicating that D r increases as F c increases, but the increased rate gradually decreases. The final D r in HF + FS and HF + TS are almost the same, thus the existence of well c has little effect on hydrate decomposition.   Figure 12 depicts the variation of Q p in HF + TS and HF + FS at different F c . In HF + TS, Q p is lower than 50 m 3 /d and gradually decreases when the F c is 5 D·cm. When F c is 10 D·cm, Q p shows a downward trend in the first 750 d as the gas released by thermal stimulation does not migrate to the well p . Subsequently, Q p increases and peaks at 85 m 3 /d. When F c is 20 and 40 D·cm, the gas migration speed is faster; thus, Q p starts to rise at 100 and 45 d, showing peak values of 143 and 232 m 3 /d, respectively. In HF + FS, Q p is higher and shows a period of increase even at a low F c value. Figure 13 Figure 14 displays the variation of V l , V c , and R e with F c . Notably, V l and V c initially increase and subsequently decrease with the increase in F c . This can be explained by the impacts of F c on hydrate decomposition and gas migration. In particular, increasing F c in the range of 5-10 D·cm significantly improved the hydrate decomposition efficiency, while the gas migration capacity was insufficient, thereby exacerbating the CH 4 leakage. When F c > 20 D·cm, the improvement of the gas migration capacity is more significant than the increase in Q r ; thus, more gas was produced by well p , which decreased the V l and V c values. However, a leakage volume of 11,236 m 3 still exists at the highest F c , indicating that the reduction of V l by increasing the F c is unsatisfactory. Additionally, V c was higher than V l , indicating that all escaped gas was recovered by well c . R e increases with increasing F c while the increased rate decreases gradually. Additionally, the R e values in HF + FS increased by 0.31, 0.21, 0.14, and 0.10 compared to those in HF + TS at F c of 5, 10, 20, and 40 D·cm, respectively. This indicated that the contribution of well c to R e decreases as F c increases.

The Effect of Well Spacing
Based on the simulation results of scenarios 04, 06-08, and 17-20#, we compared the impacts of W s on the production potential of HF + TS and HF + FS. Figure 15 displays the variation of Q r in HF + TS and HF + FS at different W s . In the first 450 days, Q r decreased as W s increased. This is attributed to the lower pressure gradient in the fracture, which reduces the fluid flow velocity. Under the same W s , the Q r in HF + FS in the first 450 days was higher than that in HF + TS, but it dropped sharply after 450 days. Figure 16 Figure 17 shows the variation of Q p in HF + TS and HF + FS at different W s . When W s is 80 m, Q p increased the fastest at the initial production stage, while it declines the most rapidly after 450 days. This is because when W s is too small, production wells are flooded earlier and there is not enough hydrate to maintain long-term production. In HF + FS, the variation in Q p presents three phases. However, in HF + TS, the pressure gradient in the fracture was insufficient for gas migration when W s ≥120 m, thus, Q p is lower and there is no rising phase. The final V p values in HF + FS were 18.62, 20.75, 21.19, and 18.70 × 10 4 m 3 at W s values of 80, 100, 120, and 140 m, respectively. These values were 1.05, 1.19, 1.47, and 1.67 times those in HF + TS (Figure 18), indicating that the enhancement of the V p increases as W s increases. Additionally, the reasonable W s should be 100-120 m because Vp under this condition is more attractive   Figure 19 depicts the variation of V l , V c , and R e with W s . Notably, V l and V c increase with increasing W s , indicating that the methane leakage is aggravated after increasing the W s . Furthermore, V c is higher than V l even at W s of 140 m, indicating no CH 4 leakage in HF + FS. The R e values in HF + TS and HF + FS are 0.91 and 1, 0.79 and 0.93, 0.6 and 0.86, and 0.54 and 0.79 at W s values of 80, 100, 120, and 140 m, respectively. Thus, with increasing W s , more gas is trapped, and the contribution of well c to R e is greater. Figure 19. The variation of V l , V c , and R e with W s .

Variance Analysis
Analysis of variance (ANOVA) can examine the significance of the difference between the means of the factors [50]. In this study, a two-factor mixed variance model is established through the "aov" function in the "R language" (https://www.r-project.org/ (accessed on 18 June 2022)). The V p , V c , and R e values of 300, 600, 900, 1200, 1500, and 1800 d in HF + FS are considered evaluation indicators, and F c and W s are the independent variables. Thus, the significance of F c and W s are obtained, as shown in Figure 20. Additionally, the main effect and interaction effect of F c and W s were visualized, as shown in Figure 21.

Gas Production
As shown in Figure 20a, the significance of F c to V p is highest in the entire production period, and the significance level of W s is lower and gradually decreases as the production progresses. Thus, F c has a more significant impact on V p . Moreover, F c has a positive effect on V p , whereas W s has a negative effect on it, except at W s of 80 m (Figure 21a-c). Thus, it is necessary to further improve F c to overcome the negative effects of W s to enhance production. For the current simulations, the average Q p per unit horizontal section length is higher than 100 m 3 /d when F c ≥ 20 D·cm, and a horizontal section length of 600 m can reach the commercial threshold of 60,000 m 3 /d [25,48,49]. Figure 20b exhibits the significance level of F c to V c is highest at 300 d, which subsequently becomes insignificant at 900 d. The significance level of W s increases gradually and reaches the maximum at 1200 d. Thus, the initial and final V c values are mainly controlled by F c and W s , respectively. As shown in Figure 21d, F c has a positive effect on V c at 300 d, indicating that the rapid decomposition of hydrate caused by the fracture increased the risk of CH 4 leakage. After 900 d (Figure 21e,f), W s exhibited a positive effect on V c . Furthermore, the effect of F c on V c changed from positive to negative as the production progresses. This is because the thermal decomposition front gradually moves away from well c , and more gas is produced when F c is higher. Therefore, a higher F c reduces V c , while the necessity of deploying well c increases under larger W s .

Recovery Efficiency
As shown in Figure 20c, both F c and W s are insignificant to R e at the initial stage, and their significance levels were enhanced as the production progresses. Moreover, the significance level of W s is higher than that of F s , indicating that W s have a more significant impact on R e . Additionally, the negative effect of W s on R e and the positive effect of F c on R e become increasingly clear as the production progresses (Figure 21g-i). Therefore, a higher R e can be obtained with a smaller W s and a greater F c , e.g., R e values were higher than 0.85 for W s ≤ 120 m and F c ≥ 20 D·cm.

Discussion
The previous study has obtained the optimal development plan of FS for nonhydraulically fractured reservoirs [29]. To evaluate the production potential of the proposed development mode in this study, the development indicators of the two are compared, as shown in Figure 22. The fracture conductivity and well spacing of the proposed development mode used for comparison here are 40 D·cm and 120 m, respectively. It can be seen that the well spacing is 50 m in a non-hydraulically fractured reservoir, which is only 41.7% of that in a hydraulically fractured reservoir. The cumulative gas production in hydraulically fractured reservoir is 243.2 × 10 3 m 3 , which is 3.1 times of non-hydraulically fractured reservoir. The volume of captured gas increases from 22.9 to 43.1 × 10 3 m 3 , indicating that there is a greater risk of gas escape in hydraulically fractured reservoir. Additionally, the recovery efficiency increases by 6.3% after fracturing, and the horizontal well section length required to reach the commercial gas production level is reduced from 1362 m to 444 m. Therefore, the proposed development mode is helpful to reduce drilling costs and improve the economics of challenging hydrate development. It is important to note that the present simulations were conducted under an ideal fracture morphology as field trials of hydraulic fracturing in hydrate reservoirs have not yet been carried out. In fact, for weakly consolidated hydrate reservoirs, it is generally considered that hydraulic fracturing is challenging, and the fracture propagation patterns are complex due to the inhomogeneous hydrates in sediments [51][52][53]. Thus, it is necessary to further evaluate the effect of fracturing morphology on reservoir stimulation and methane leakage. Additionally, various injection-production well patterns, such as two-, three-, five-, seven-, and nine-spot, have been designed in conventional oil and gas development, and there is an optimal matching relationship between fracture distribution and injection-production system [54][55][56]. Unlike conventional oil and gas, the design of hydrate development mode needs to consider the methane leakage due to the lack of impermeable boundaries. Furthermore, the methane escape behavior is affected by reservoir-cover conditions and injection-production parameters. Consequently, to safely and efficiently recover methane from challenging hydrates, the applicability of different development modes to challenging hydrates with various occurrence conditions should be further explored.

Conclusions
In this study, we compared the production potential of the three-spot and five-spot well patterns for hydraulically fractured CH-kb reservoirs. The effects of fracture conductivity and well spacing on methane leakage and gas production were investigated. The main conclusions are as follows.

1.
For the three-spot well pattern, there was considerable methane leakage into the permeable boundary layers. The leakage volume increases as the well spacing increases and reaches 46,341 m 3 when the well spacing is 140 m. The impact of fracture conductivity on methane leakage shows two opposite stages: when the fracture conductivity is lower than 20 D·cm, the leakage volume increases as the fracture conductivity increases; when the fracture conductivity is higher than 20 D·cm, the opposite is true. However, enhancing the fracture conductivity can not prevent methane leakage, as a leakage volume of 11,236 m 3 was observed at the greatest fracture conductivity of 40 D·cm.

2.
The proposed five-spot well pattern combined with hydraulic fracturing could realize the safe development of CH-kb under a large well spacing. Even at the well spacing of 140 m, there is still no methane leakage. Additionally, increasing the well spacing increased the risk of methane leakage. Consequently, deploying capture wells is crucial for expanding safe well spacing.

3.
Compared with the three-spot, the five-spot further increases the gas production and recovery and is expected to meet the commercial production threshold of hydrates. In the current simulations, the five-spot maximizes total gas production and recovery efficiency by 87 and 31%, respectively, compared to the three-spot. The average gas production rate reached 142 m 3 /d per unit horizontal section length, and a horizontal section length of 425 m was estimated to meet the commercial threshold.

4.
The development mode of hydraulic fracturing combined with the five-spot well pattern enhances the production potential of CH-kb while addressing methane leakage, which is safe and efficient for the development of low-permeable hydrate reservoirs with permeable boundary layers. In particular, a comprehensive analysis of gas production, methane capture, and gas recovery indicates that the CH-kb hydrate reservoir at the SH2 site in the Shenhu area is suitable for a well spacing of 100-120 m and fracture conductivity greater than 20 D·cm. Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.