Numerical Evaluation of Commingled Production Potential of Marine Multilayered Gas Hydrate Reservoirs Using Fractured Horizontal Wells and Thermal Fluid Injection

Abstract

Multilayered reservoirs with coexisting free gas and hydrates are primary targets for commercialization, nevertheless, the extremely low permeability greatly limits their extraction efficiency. Herein, multilayer commingled production using horizontal wells stimulated by hydraulic fracturing and thermal fluid injection was proposed to enhance productivity, and the effects of key factors on co-production performance were numerically examined, with the reservoir located in the Shenhu Area as the geological background. The results indicated that due to severe interlayer contradictions, the stimulation capabilities of using fracturing or thermal fluid injection alone were limited, in particular, the extraction of hydrates severely lagged behind. However, their combination exhibited tantalizing productivity due to strengthened inter-well interaction. Reducing the fracture spacing was more effective than increasing fracture conductivity in shortening the production cycle, and intensive fractures with adequate flow capacity were suggested for gas enhancement and water control. When the fracture spacing was reduced from 30 to 5 m and the fracture conductivity increased from 10 to 100 D·cm, the horizontal section length for commercial production (average daily gas production of 50,000 m³ and recovery ratio of 0.7) was reduced from 1758 to 146 m, which is lower than the on-site horizontal section length of 250–300 m. Therefore, the proposed development mode is promising for the commingled production of gas and hydrates.

1. Introduction

Natural gas hydrate (NGH) is a cage-like crystallized compound consisting of gas (mainly CH₄) and water that is widespread in submarine sediments and permafrost regions [1]. The organic carbon content in naturally occurring NGH is estimated to be about 10⁴ billion tons, approximately one-third of the world’s total. As a result, NGH is widely recognized as the ideal new clean energy resource for the 21st century [2].

Gas recovery from NGH faces enormous challenges due to the harsh environment, complex phase behavior, and potential geo-environmental risks. Currently, only 10 production trials have been carried out at the Mallik site in Canada [3], the North Slope of Alaska [4], the Nankai Trough [5], the Qilian Mountains [6], and the Shenhu Area [7]. The production trials in the Shenhu Area especially demonstrated the feasibility of gas extraction from oceanic clayey-silty sediments. This challenging storage environment contains more than 90% of the world’s NGH resources and is generally considered the most difficult to extract due to its low effective permeability [8]. However, since there are three sublayers, that is, the hydrate-bearing layer (HBL), three-phase layer (TPL), and free gas layer (FGL), and the use of horizontal well technology, the productivity is the highest. Additionally, such multilayered reservoirs with coexisting NGH and free gas are quite prevalent and have been found at Blake Ridge [9], Green Canyon [10], and Offshore Oregon [11], and are, therefore, generally recognized as the “sweet spots” and primary targets for commercialization.

Recently, aiming for the commercialization of gas hydrates, various development modes have been raised, including the following: i. thermal stimulation to provide sufficient energy for hydrate dissociation, with heating methods including massless injection in the form of electrical or microwave heating [12,13], as well as thermal fluid injection (mainly hot water) in the form of single-well huff and puff [14] and injection-production well patterns [15]; ii. using horizontal wells to increase the drainage area, such as radially branched horizontal wells [16]; iii. implementing reservoir reconstruction such as enlarged permeable well walls [17] and hydraulic fracturing [18] to improve reservoir permeability; iv. injecting gases (mainly CO₂) that are more likely to form hydrates to replace CH₄ molecules in NGH cages [19]; v. injecting chemical inhibitors like brine to promote NGH dissociation [20]; vi. integrating multiple methods, such as combining fracturing and multi-well systems [21], and using thermal stimulation and chemical inhibitors [22]. Nevertheless, these studies mainly focus on the production dynamics of reservoirs with a single HBL, with few simulations exploring the commingled production of multilayers. Yu et al. (2021) found that the gas contribution of each sublayer was FGL, HBL, and TPL in descending order for a vertical well [23]. Subsequently, the exploitation dynamics of the horizontal well with multiple artificial fractures were further simulated, and the results indicated that horizontal fractures contributed more to productivity [24]. Ma et al. (2021) examined the influences of well location and found that deploying the horizontal section in the HBL was beneficial for long-term production [25]. Yin et al. (2022) suggested that the low permeability of the HBL inhibits pressure and temperature transfer, leading to inefficient hydrate decomposition [26]. Overall, the above studies preliminarily proved the feasibility of implementing multilayer co-production using horizontal wells. However, several crucial issues need further research: firstly, the capabilities of various means of stimulation like hot water injection and fracturing are poorly understood; secondly, the response mechanism of co-production dynamics to various stimulation factors is still unclear; and lastly, the commercial development potential of multilayer co-production has not yet been effectively assessed.

This study aims to investigate the above issues through simulations of commingled production in the multilayered reservoir located in the Shenhu Area. Firstly, co-production performance using horizontal wells stimulated by hydraulic fracturing and/or thermal fluid injection was compared, followed by the effects of fractured sublayers, fracture spacing, and fracture conductivity on co-production dynamics, including hydrate decomposition and gas–water production; and finally, through variance analysis and commercial potential evaluation, the optimal multilayer co-production plan was obtained. The findings of this study help to understand the peculiar multilayer co-production behavior and stimulation mechanism, as well as provide references for the implementation of reservoir stimulation and the design of co-production schemes. In addition, the proposed multilayer co-production with the combined stimulation of horizontal well fracturing and hot water injection offers an attractive development mode for this kind of low-permeability multilayered “sweet-spot” gas hydrate reservoir, which is conducive to promoting the industrialization process of marine-challenging hydrates.

2. Numerical Modeling

2.1. Geoengineering Background

The Shenhu Area, located in the South China Sea, is the main NGH-rich area in China. Since 2007, three gas hydrate geological surveys have been carried out, during which samples with a hydrate saturation of 20–48% and CH₄ content of more than 99.7% were recovered [27]. Recent research on the deposition system has demonstrated the coexistence of hydrate and free gas by analyzing geology and logging data [28].

2.2. Development Mode Design

Here, the multilayered reservoir with a seawater depth of 1225.23 m at the test site was taken as the research object, as shown in Figure 1. The target production layers are located at 207.8–297 mbsf. The thicknesses of the overlayer (OL) and underlayer (UL) were set to 30 m to accurately simulate the reservoir temperature and pressure evolution [29]. A three-spot well pattern was used with all horizontal sections extending along the Y-axis. Two depressurization wells were arranged at the top of the TPL for the co-production of free gas and NGH dissociation gas, while the middle injection well was deployed at the top of the HBL to maximize the hot water sweep area. Reservoir properties were assumed to be homogeneous in the horizontal since there is no available data for reference, and thus, only a 30 m horizontal well section was considered for simulation.

It was assumed that multiple vertical fractures were formed by segmental fracturing, as shown in Figure 2. The well spacing was set to 30 m for inter-well connectivity. To reveal the impact of fractured sublayers on co-production performance, three cases were considered, that is, fractured HBL (Figure 2a), fractured HBL and TPL (Figure 2b), and fractured HBL, TPL, and FGL (Figure 2c). In addition, the value of fracture spacing was set to four levels: 30, 15, 10, and 5 m, and thus, a 30 m horizontal section contained two, three, four, and seven fractures, respectively, as shown in Figure 2d–g.

2.3. Numerical Simulator

Several numerical simulators have been developed for NGH exploitation, such as CMG-STARS, Tough + Hydrate (T + H), and HydrResSim. Among them, T + H is preferred because of its powerful functions in mine scale simulation. Moreover, its accuracy has been fully examined [30] and thus was employed here.

NGH exploitation is accompanied by phase transition, as follows

$${\mathrm{C}\mathrm{H}}_4·5.75{\mathrm{H}}_2\mathrm{O}={\mathrm{C}\mathrm{H}}_4+5.75{\mathrm{H}}_2\mathrm{O}$$

(1)

The mass and energy balances in each domain are expressed as [31]

$$\frac{d}{dt}{\int }_{{\mathrm{V}}_{\mathrm{n}}}^{ }{M}^{\kappa }d\mathrm{V}={\int }_{{\mathsf{\Gamma }}_{\mathrm{n}}}^{ }{F}^{\kappa }dV·\mathrm{n}d\mathsf{\Gamma }+{\int }_{{\mathrm{V}}_{\mathrm{n}}}^{ }{q}^{\kappa }d\mathrm{V},\kappa \equiv m\left({\mathrm{C}\mathrm{H}}_4\right),w(\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}),i(\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{i}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{o}\mathrm{r}),\theta (\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t})$$

(2)

The mass accumulation term is given by

$$M^{\kappa }=\sum _{\beta =H,A,G,I}\varphi {\rho }_{\beta }{S}_{\beta }{X}_{\beta }^{\kappa }, \beta \equiv H\left(\mathrm{h}\mathrm{y}\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\right), A\left(\mathrm{a}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}\right), G\left(\mathrm{g}\mathrm{a}\mathrm{s}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}\right), I(\mathrm{i}\mathrm{c}\mathrm{e})$$

(3)

The energy accumulation term consists of the matrix energy, energies of all phases, and NGH dissociation activation energy, as follows

$$M^{\theta }=\left(1-\varphi \right){\rho }_{sed}{C}_{sed}\mathrm{T}+\sum _{\beta =G,A,H,I}\varphi {\rho }_{\beta }{S}_{\beta }{U}_{\beta }+E_{dis}$$

(4)

In the code, Darcy’s law is employed, while taking into account the diffusion effect of the gaseous and multiphase fluid flow, and is described by

$$F_A^{\kappa }=-k\frac{k_{rA}{\rho }_A{X}_A^{\kappa }}{{\mu }_A}\left(\nabla P_A-{\rho }_{A}g\right)$$

(5)

$$F_G^{\kappa }=-k\frac{k_{rG}{\rho }_G{X}_G^{\kappa }}{{\mu }_G}\left(1+\frac{\omega }{P_G}\right)\left(\nabla P_G-{\rho }_{G}g\right)$$

(6)

The relative permeabilities of the gaseous and aqueous are given by [32]

$$k_{rA}={\left(\frac{S_A-S_{irA}}{1-S_{irA}}\right)}^{n_A}$$

(7)

$$k_{rG}={\left(\frac{S_G-S_{irG}}{1-S_{irA}}\right)}^{n_G}$$

(8)

The aqueous pressure considering capillary pressure is expressed as

$$P_A=P_G+P_{cap}$$

(9)

The capillary pressure in the gas hydrate system is described by [33]

$$P_{cap}=-P_{thr}{\left[{\left(\frac{S_A-S_{irA}}{1-S_{irA}}\right)}^{-\frac{1}{\lambda }}-1\right]}^{1-\lambda }$$

(10)

The energy flux is expressed as

$$F^{\theta }=f_{\sigma }{\sigma }_0\nabla {\mathrm{T}}^4-\left[\left(1-\varphi K_{sed}+\varphi \sum _{\beta =G,A,H,I}{S}_{\beta }{K}_{\beta }\right)\right]\nabla \mathrm{T}+\sum _{\beta =G,A}{h}_{\beta }{F}_{\beta }$$

(11)

The mass and energy source-sink terms are expressed as follows

$$q^{\kappa }=\sum _{\beta =G,A}{X}_{\beta }^{\kappa }{q}_{\beta }$$

(12)

$$q^{\theta }=q_d+\sum _{\beta =G,A}{h}_{\beta }{q}_{\beta }$$

(13)

2.4. Initial and Boundary Conditions

The pore seawater was assumed to be connected because of the open boundaries, and thus, the pore pressure gradient was calculated by

$$P_{por}=P_{atm}+{\rho }_{sw}g\left(h_0+z\right)$$

(14)

In addition, the temperature field was initialized based on the following equation

$$\mathrm{T}={\mathrm{T}}_{\mathrm{s}\mathrm{f}}+∆\mathrm{T}·z$$

(15)

Dirichlet boundary conditions were imposed at upper and lower boundaries, with constant pressures and temperatures there. Additionally, the porosity and intrinsic permeability of the OL and UL were not available, and are, therefore, presumed to be the same as those of the adjacent sublayers. Notably, the temperature–pressure transmission and fluid flow induced by OL and UL are greatly weakened because of the use of hydraulic fracturing and thermal fluid injection, making the impact of this simplification on the simulation results very limited. The model parameters are presented in

Table 1. Model parameters.

Parameters	Value & Unit
Porosity of OL, HBL, TPL, FGL, and UL	0.373, 0.373, 0.346, 0.347, and 0.347
Intrinsic permeability of OL, HBL, TPL, FGL, and UL	2.38, 2.38, 6.63, 6.8, and 6.8 mD
Mean hydrate saturation in HBL and TPL	31% and 11.7%
Mean gas saturation in TPL and FGL	13.2% and 9.3%
Gas composition	CH4 100%
Seawater salinity	0.03
Dry and wet thermal conductivity (all deposits)	1 and 3.1 W/m/K
Sediment density	2650 kg/m3
Sediment specific heat	1000 J/K/kg
Production and injection pressure	3 and 16 MPa
Hot water temperature	60 °C
${\rho }_{sw}$	1024 kg/m3
${\mathrm{T}}_{\mathrm{s}\mathrm{f}}$	3.6 °C
$∆\mathrm{T}$	46 °C/km
$P_{thr}$	1 × 105 Pa
$\lambda$	0.45
$n_A$	3.50
$n_G$	3.50
$S_{irA}$	0.30
$S_{irG}$	0.03

.

2.5. Model Validation

A simulation was conducted on the unit-length horizontal section for model validation, and the simulation data are shown in Figure 3. As can be seen, 3312 m³ of gas was recovered from the unit-length horizontal section over 30 days. This results in 861,400 m³ of gas for a horizontal section length of 260 m, within the range of the on-site horizontal section length of 250–300 m [34], proving the accuracy of the model.

2.6. Simulation Scenarios

The co-production performances of four development modes were compared, that is, depressurization without stimulation (D), stimulation using fracturing (D+F) or thermal fluid injection (D+T), or both (D+T+F). Additionally, referring to the existing research, the fracture conductivity was set to 10–100 D·cm [21]. Ultimately, the factors examined include stimulation methods, fractured sublayers, fracture spacing, and fracture conductivity, with simulation scenarios presented in

Table 2. Simulation scenarios.

Scenario No.	Development Mode	Fractured Sublayers	${\mathbf{C}}_{\mathbf{F}}$ (D·cm)	${\mathbf{S}}_{\mathbf{F}}$ (m)	Scenario No.	Development Mode	Fractured Sublayers	${\mathbf{C}}_{\mathbf{F}}$ (D·cm)	${\mathbf{S}}_{\mathbf{F}}$ (m)
01#	D	/	/	/	12#	D+T+F	HBL, TPL	25	15
02#	D+T	/	/	/	13#	D+T+F	HBL, TPL	25	10
03#	D+F	HBL	10	30	14#	D+T+F	HBL, TPL	25	5
04#	D+T+F	HBL	10	30	15#	D+T+F	HBL, TPL	50	30
05#	D+T+F	HBL, TPL	10	30	16#	D+T+F	HBL, TPL	50	15
06#	D+T+F	HBL, TPL, FGL	10	30	17#	D+T+F	HBL, TPL	50	10
07#	D+T+F	HBL, TPL	10	30	18#	D+T+F	HBL, TPL	50	5
08#	D+T+F	HBL, TPL	10	15	19#	D+T+F	HBL, TPL	100	30
09#	D+T+F	HBL, TPL	10	10	20#	D+T+F	HBL, TPL	100	15
10#	D+T+F	HBL, TPL	10	5	21#	D+T+F	HBL, TPL	100	10
11#	D+T+F	HBL, TPL	25	30	22#	D+T+F	HBL, TPL	100	5

.

3. Results and Discussion

3.1. Comparison of Different Development Modes

The co-production performance under different development modes was compared based on scenario numbers 01~04#.

3.1.1. NGH Dissociation

Figure 4 depicts the evolution of ${\mathrm{Q}}_{\mathrm{R}}$ and ${\mathrm{R}}_{\mathrm{D}}$ under different development modes. The ${\mathrm{Q}}_{\mathrm{R}}$ in the D mode exhibited an initial sharp drop and a subsequent (~300 days) slow decrease. The mechanism for this is twofold: first, the NGH dissociation front gradually moved away from the wellbores (Figure 5a–c); and second, the reservoir temperature decreased significantly (Figure 6a–c) and cannot provide enough heat for hydrate dissociation. In addition, the hydrate around the middle horizontal well was barely decomposed because boundary water intrusion hinders pressure drop transmission [35,36]. Owing to the greater permeability and the critical steady state of the hydrates, a larger hydrate decomposition radius was observed in the TPL (Figure 5c). In the D+T mode, the improvement in ${\mathrm{Q}}_{\mathrm{R}}$ was not significant as hot water injection and migration were suppressed by the HBL’s low permeability (Figure 5d–f), and substantial hot water leakage into the OL can be observed in Figure 6d–f, resulting in only a 15% increase in the final ${\mathrm{R}}_{\mathrm{D}}$ (Figure 4). In the D+F mode, the maximum ${\mathrm{Q}}_{\mathrm{R}}$ was increased to 3974 m³/d as rapid NGH dissociation occurred along fractures (Figure 5g–i), and the reservoir temperature dropped significantly due to the limited sensible heat (Figure 6g–i), resulting in a sharp decrease in ${\mathrm{Q}}_{\mathrm{R}}$ [37]. At 1200 days, the hydrates in the TPL were almost completely decomposed, whereas a substantial amount of hydrates was still reserved in the HBL (Figure 5i). Thus, the extraction efficiency of the HBL was still severely lagging even after fracturing. In the D+T+F mode, ${\mathrm{Q}}_{\mathrm{R}}$ was greatly improved with a high ${\mathrm{Q}}_{\mathrm{R}}$ (>1000 m³/d) phase (~1000 days). This is because the fractures intensify the hot water injection and provide high-conductive paths for fluid transport [38,39], thereby enlarging the hot water sweep area (Figure 6j–l) and facilitating hydrate decomposition (Figure 5j–l). The final ${\mathrm{R}}_{\mathrm{D}}$ reached 99.5%, 3.98, 2.49, and 2.24 times that in D, D+T, and D+F, respectively. Therefore, hydrate decomposition efficiency is greatly improved by combining thermal stimulation and fracturing.

3.1.2. Gas and Water Production

Figure 7 depicts the evolution of ${\mathrm{Q}}_{\mathrm{G}}$ and ${\mathrm{V}}_{\mathrm{G}}$ under different development modes. ${\mathrm{Q}}_{\mathrm{G}}$ and ${\mathrm{V}}_{\mathrm{G}}$ in the D and D+T modes were almost the same because neither thermal stimulation nor depressurization at the top of the HBL was ineffective. The main gas source layers were the TPL and FGL, with very limited contribution from the HBL (Figure 8a–f). In the D+F mode, ${\mathrm{Q}}_{\mathrm{G}}$ was increased as the decomposition gas was produced (Figure 8g–i), but the stimulation capacity of the fractures was still limited. Unlike the decay of ${\mathrm{Q}}_{\mathrm{G}}$ in the other three scenarios, there was a rise phase (~280 days) after the initial sharp drop of ${\mathrm{Q}}_{\mathrm{G}}$ (Figure 7) due to the contribution of thermal decomposition gas (Figure 8j–l). When the hot water breaks through into the production wells along the fractures (~800 days), ${\mathrm{Q}}_{\mathrm{G}}$ drops dramatically because the hot water in the fracture impedes gas output. The final ${\mathrm{V}}_{\mathrm{G}}$ reached 214 × 10⁴ m³, 2.15, 2.22, and 1.79 times that in the D, D+T, and D+F modes, respectively, suggesting that the D+T+F mode is the most promising.

Figure 9 depicts the evolution of ${\mathrm{V}}_{\mathrm{W}}$ and ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ under different development modes. The initial ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ in the D and D+T modes was higher, but afterward, it dropped sharply because of the decreases in ${\mathrm{Q}}_{\mathrm{G}}$ and the increase in ${\mathrm{V}}_{\mathrm{W}}$. In the D+F mode, the initial ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ was the lowest, while there was an initial increase in ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ due to the stimulation of fractures. The ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ decay rate in the D+T+F mode was the lowest due to the relatively stable ${\mathrm{Q}}_{\mathrm{G}}$. Overall, when ${\mathrm{V}}_{\mathrm{G}}$ was 130 × 10⁴ m³, the ${\mathrm{V}}_{\mathrm{W}}$ in the D+T+F mode was the lowest, and the corresponding ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ reached 15.72, which was 5.40, 2.69, and 2.61 times that in the D, D+T, and D+F modes, respectively. Consequently, it also manifests the desirable water production performance in the D+T+F mode.

3.2. Effect of Fractured Sublayers

The effect of fractured sublayers was analyzed based on scenario numbers 04~06#.

3.2.1. NGH Dissociation

Figure 10 depicts the evolution of ${\mathrm{Q}}_{\mathrm{R}}$ and ${\mathrm{R}}_{\mathrm{D}}$ under various fractured sublayers. The ${\mathrm{Q}}_{\mathrm{R}}$ in the first 500 days was improved by fracturing the TPL as hydrate dissociation was encouraged in this sublayer (Figure 11a,b), whereas fracturing the FGL hardly affects the ${\mathrm{Q}}_{\mathrm{R}}$ and ${\mathrm{R}}_{\mathrm{D}}$ (Figure 11b,c). Consequently, simultaneous fracturing of the HBL and TPL is a good plan, which can completely decompose the hydrate within 1767 days.

3.2.2. Gas and Water Production

Figure 12 and Figure 13 depict the evolution of ${\mathrm{Q}}_{\mathrm{G}}$, ${\mathrm{V}}_{\mathrm{G}}$, ${\mathrm{V}}_{\mathrm{W}}$, and ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ under different fractured sublayers. The initial ${\mathrm{Q}}_{\mathrm{G}}$ increased from 2565 to 3278 m³/d by fracturing the TPL due to the increased ${\mathrm{Q}}_{\mathrm{R}}$ and the contribution of free gas (Figure 11d,e). Whereas fracturing the FGL reduced ${\mathrm{Q}}_{\mathrm{G}}$ because bottom-water intrusion decreases the gas transport capacity of the fractures (Figure 11e,f). This can also be confirmed by Figure 13, where the final ${\mathrm{V}}_{\mathrm{W}}$ increased from 18 to 21 × 10⁴ m³ and the final ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ decreased from 10.94 to 9.45 after fracturing the FGL. In addition, the initial ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ was increased after fracturing the TPL and the final ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ was acceptable. As a whole, fracturing the HBL and TPL is the optimal scheme.

3.3. Effect of Fracture Conductivity and Fracture Spacing

${\mathrm{C}}_{\mathrm{F}}$ reflects the rate of fluid transport within fractures, while ${\mathrm{S}}_{\mathrm{F}}$ affects the distribution and effective surface area of fractures, and both of them influence fluid migration and production. Therefore, a comparison of the effects of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ on co-production dynamics is required to obtain a desirable stimulation scheme. Here, based on scenario numbers 07~11#, 15#, and 19#, the effects of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ were analyzed.

3.3.1. NGH Dissociation

Figure 14 depicts the evolution of ${\mathrm{Q}}_{\mathrm{R}}$ and ${\mathrm{R}}_{\mathrm{D}}$ at different ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ values. The initial ${\mathrm{Q}}_{\mathrm{R}}$ increased from 5132 to 33,797 m³/d as ${\mathrm{C}}_{\mathrm{F}}$ increased from 10 to 100 D·cm, whereas it was 13,765 m³/d at an ${\mathrm{S}}_{\mathrm{F}}$ of 5 m, indicating that increasing ${\mathrm{C}}_{\mathrm{F}}$ was more effective. A higher decay rate of ${\mathrm{Q}}_{\mathrm{R}}$ was observed at a higher ${\mathrm{C}}_{\mathrm{F}}$ because of the premature thermal breakthrough [38]. In addition, the enhancement of ${\mathrm{Q}}_{\mathrm{R}}$ was not significant (Figure 15a–d) when ${\mathrm{C}}_{\mathrm{F}}$ was higher than 25 D·cm because of insufficient heat exchange. Different from increasing ${\mathrm{C}}_{\mathrm{F}}$, the reduction of ${\mathrm{S}}_{\mathrm{F}}$ enhanced heating efficiency due to the increased hot water sweep area, and thus exhibited a higher ${\mathrm{Q}}_{\mathrm{R}}$ (Figure 15e–g). The hydrate dissociation cycle was reduced from 1768 to 1220 days as ${\mathrm{C}}_{\mathrm{F}}$ increased from 10 to 100 D·cm, whereas it was 645 days at an ${\mathrm{S}}_{\mathrm{F}}$ of 5 m. Consequently, compared with increasing ${\mathrm{C}}_{\mathrm{F}}$, the decrease in ${\mathrm{S}}_{\mathrm{F}}$ is more effective for hydrate dissociation.

3.3.2. Gas and Water Production

Figure 16 and Figure 17 depict the evolution of ${\mathrm{Q}}_{\mathrm{G}}$, ${\mathrm{V}}_{\mathrm{G}}$, ${\mathrm{V}}_{\mathrm{W}}$, and ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ at different ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ values. Increasing ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ mainly increased ${\mathrm{Q}}_{\mathrm{G}}$ in the first 400 days. The difference is that ${\mathrm{Q}}_{\mathrm{G}}$ decreased sharply at higher ${\mathrm{C}}_{\mathrm{F}}$ values, while it was relatively stable at different ${\mathrm{S}}_{\mathrm{F}}$ values. Moreover, ${\mathrm{V}}_{\mathrm{W}}$ increased significantly as ${\mathrm{C}}_{\mathrm{F}}$ increased, while it was barely affected by ${\mathrm{S}}_{\mathrm{F}}$. Overall, the time and ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ for a ${\mathrm{V}}_{\mathrm{G}}$ of 200 × 10⁴ m³ was reduced to 550 days and 3.24 as ${\mathrm{C}}_{\mathrm{F}}$ increased to 100 D·cm, whereas they were 535 days and 12.26 at an ${\mathrm{S}}_{\mathrm{F}}$ of 5 m. As a result, both increasing ${\mathrm{C}}_{\mathrm{F}}$ and decreasing ${\mathrm{S}}_{\mathrm{F}}$ help gas recovery, with no increase in water production even at a low ${\mathrm{S}}_{\mathrm{F}}$ value.

3.4. Vanrican Analysis

Analysis of variance was employed based on the “aov” function in the “R language” (https://www.r-project.org (accessed on 18 June 2022)) to reveal the significance of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ and their interactions on co-production performance, and the significance of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ to $\overline{{\mathrm{Q}}_{\mathrm{R}}}$, ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$, and $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ are shown in Figure 18. In addition, the main effect and interaction effect were visualized using the “interractionwt” function, as shown in Figure 19.

3.4.1. Significance of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ to $\overline{{\mathrm{Q}}_{\mathrm{R}}}$

The p-values of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ in relation to $\overline{{\mathrm{Q}}_{\mathrm{R}}}$ increase as ${\mathrm{R}}_{\mathrm{D}}$ increases (Figure 18a), indicating their significance gradually decreases as production progresses. When ${\mathrm{R}}_{\mathrm{D}}$ was 0.1, ${\mathrm{C}}_{\mathrm{F}}$ was more significant, whereas the significance of ${\mathrm{S}}_{\mathrm{F}}$ was higher when ${\mathrm{R}}_{\mathrm{D}}$ reached 0.5. This indicates that the initial ${\mathrm{R}}_{\mathrm{D}}$ is primarily controlled by ${\mathrm{C}}_{\mathrm{F}}$, as it directly impacts the rate of fluids (such as gas) flowing through the fractures, thereby greatly affecting the initial gas production rate; and for ${\mathrm{S}}_{\mathrm{F}}$, its impact on the final ${\mathrm{R}}_{\mathrm{D}}$ is greater because reducing fracture spacing allows for more efficient drainage of the remaining gas, leading to a more significant impact on the final level of reservoir depletion. In addition, ${\mathrm{C}}_{\mathrm{F}}$ has a positive effect on $\overline{{\mathrm{Q}}_{\mathrm{R}}}$ due to the increased hot water flow rate at a lower ${\mathrm{C}}_{\mathrm{F}}$, while ${\mathrm{S}}_{\mathrm{F}}$ has a negative effect (Figure 19a–c) due to the increased thermal fluid sweep area at a smaller ${\mathrm{S}}_{\mathrm{F}}$. Therefore, NGH dissociation requires a high ${\mathrm{C}}_{\mathrm{F}}$ and a lower ${\mathrm{S}}_{\mathrm{F}}$.

3.4.2. Significance of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ to $\overline{{\mathrm{Q}}_{\mathrm{G}}}$

The p-values of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ in relation to $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ first increase and then decrease as ${\mathrm{R}}_{\mathrm{e}}$ increases (Figure 18b). A maximum p-value difference occurred at an ${\mathrm{R}}_{\mathrm{e}}$ of 0.1 because the initial gas production was mainly controlled by ${\mathrm{C}}_{\mathrm{F}}$. When ${\mathrm{R}}_{\mathrm{e}}$ reached 0.4, the significance codes of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ were the same, suggesting that both ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ determined the subsequent gas production. The main and interaction effects of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ on $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ and $\overline{{\mathrm{Q}}_{\mathrm{R}}}$ were similar (Figure 19a–f). The difference was that their significance to $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ was greater as fractures were the primary channels for gas output. Thus, high-density, high-conductivity fractures are required to enhance co-production efficiency.

3.4.3. Significance of ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ to ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$

${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ were significant to ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ at different ${\mathrm{R}}_{\mathrm{e}}$ values (Figure 18c) because both ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ have a significant impact on gas and water production. When R_e was 0.3, ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ was the highest at a ${\mathrm{C}}_{\mathrm{F}}$ of 25 D·cm and an ${\mathrm{S}}_{\mathrm{F}}$ of 10 m (Figure 19g). Whereas ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ decreased as ${\mathrm{C}}_{\mathrm{F}}$ and ${\mathrm{S}}_{\mathrm{F}}$ increased at an ${\mathrm{R}}_{\mathrm{e}}$ of 0.5, and thus the highest ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ occurred at a ${\mathrm{C}}_{\mathrm{F}}$ of 10 D·cm and an ${\mathrm{S}}_{\mathrm{F}}$ of 5 m (Figure 19h). This is because reducing ${\mathrm{S}}_{\mathrm{F}}$ provides more gas production channels, while lowing ${\mathrm{C}}_{\mathrm{F}}$ reduces the hot water flow rate, thereby increasing gas production and reducing water production. When ${\mathrm{R}}_{\mathrm{e}}$ was 0.7, the ideal ${\mathrm{R}}_{\mathrm{G}\mathrm{W}}$ can be obtained (Figure 19i) at a ${\mathrm{C}}_{\mathrm{F}}$ of 10 D·cm and an ${\mathrm{S}}_{\mathrm{F}}$ of 5–15 m. In conclusion, if both water control and gas enhancement are considered, one option is to increase the number of fractures but at the same time control the flow capacity of the fractures.

3.5. Commercial Potential Evaluation

Commercial potential evaluation is crucial for the assessment of NGH resources and the design of development plans. According to the experience of two production tests and in conjunction with the economic evaluation, Japan suggests that the $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ of a single well for the industrialized development of marine NGH at least reach 50,000 m³ [40]. Given further integration with the receptacle ${\mathrm{R}}_{\mathrm{e}}$ value of 0.7 for water-flood development in Chinese offshore gas reservoirs [41], the commercial potential of the proposed development mode was evaluated. Specifically, when the ${\mathrm{R}}_{\mathrm{e}}$ value reaches 0.7, commercial development is possible if the $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ of a single well is higher than 50,000 m³; otherwise, it is not.

Figure 20 shows the ${\mathrm{L}}_{\mathrm{W}}$ for reaching ${\mathrm{R}}_{\mathrm{e}}$ = 0.7 and $\overline{{\mathrm{Q}}_{\mathrm{G}}}$ = 50,000 m³. The maximum ${\mathrm{L}}_{\mathrm{W}}$ is 1758 m occurring at a ${\mathrm{C}}_{\mathrm{F}}$ of 10 D·cm and an ${\mathrm{S}}_{\mathrm{F}}$ of 30 m. With the increase of ${\mathrm{C}}_{\mathrm{F}}$ to 100 D·cm and the decrease of ${\mathrm{S}}_{\mathrm{F}}$ to 5 m, it gradually decreases to 146 m. According to the on-site ${\mathrm{L}}_{\mathrm{W}}$ of 250 to 300 m, commercial development of Shenhu hydrates can be realized when ${\mathrm{C}}_{\mathrm{F}}$ is 50 D·cm and ${\mathrm{S}}_{\mathrm{F}}$ = 5 m or when ${\mathrm{C}}_{\mathrm{F}}$ = 100 D·cm and ${\mathrm{S}}_{\mathrm{F}}$ ≤ 10 m.

4. Conclusions

The co-production potential of multilayered gas hydrate reservoirs using horizontal wells stimulated by hydraulic fracturing and thermal fluid injection was numerically evaluated, and the influence mechanisms of various factors on co-production behavior were revealed. The following conclusions can be drawn:

1. The stimulation capability of fracturing or thermal fluid injection alone is limited due to severe interlayer contradictions. However, their combination exhibits tantalizing gas production as the inter-well interaction is greatly strengthened by high-conductivity fractures, including promoting pressure drop transmission, enlarging the thermal fluid sweep area, and facilitating fluid injection and gas recovery.

2. Fracturing the HBL and TPL is conducive to hydrate dissociation, while fracturing the FGL induces severe bottom-water intrusion, which greatly increases water production. As a result, fracturing the HBL and TPL simultaneously is recommended.

3. By increasing fracture conductivity and reducing fracture spacing, the hydrate dissociation cycle can be greatly shortened. Moreover, reducing fracture spacing is more effective due to the increased thermal fluid sweep area. As a result, intensive fractures with adequate flow capacity are effective for water control and gas enhancement.

4. The horizontal section length for reaching the commercial development potential is reduced from 1758 to 146 m as the fracture conductivity improved from 10 to 100 D·cm and the fracture spacing reduced from 30 to 5 m, which greatly reduces the exploitation difficulty. Therefore, fractured horizontal wells combined with thermal fluid injection are promising for these multilayered gas hydrate reservoirs.

These findings can provide valuable references for the design of co-production plans and advance the commercialization process of this multilayered gas hydrate reservoir. However, the current simulations were performed under relatively ideal geoengineering conditions due to the lack of sufficient reference data. The inhomogeneous distribution of hydrates, reservoir heterogeneity, and fracture morphology may affect thermal fluid migration, hydrate dissociation, and gas–water flow, resulting in potential biases. Therefore, further research should be conducted on these issues.