Production Simulation of Stimulated Reservoir Volume in Gas Hydrate Formation with Three-Dimensional Embedded Discrete Fracture Model

Abstract

Natural gas hydrates (NGHs) in the Shenhu area of the South China Sea are deposited in low-permeability clayey silt sediments. As a renewable energy source with such a low carbon emission, the exploitation and recovery rate of NGH make it difficult to meet industrial requirements using existing development strategies. Research into an economically rewarding method of gas hydrate development is important for sustainable energy development. Hydraulic fracturing is an effective stimulation technique to improve the fluid conductivity. In this paper, an efficient three-dimensional embedded discrete fracture model is developed to investigate the production simulation of hydraulically fractured gas hydrate reservoirs considering the stimulated reservoir volume (SRV). The proposed model is applied to a hydraulically fractured production evaluation of vertical wells, horizontal wells, and complex structural wells. To verify the feasibility of the method, three test cases are established for different well types as well as different fractures. The effects of fracture position, fracture conductivity, fracture half-length, and stimulated reservoir volume size on gas production are presented. The results show that the production enhancement in multi-stage fractured horizontal wells is obvious compared to that of vertical wells, while spiral multilateral wells are less sensitive to fractures due to the distribution of wellbore branches and perforation points. Appropriate stimulated reservoir volume size can obtain high gas production and production efficiency.

1. Introduction

As a potential clean energy source, natural gas hydrate (NGH) has been predicted to have twice as many resources distributed globally as conventional fossil fuels combined [1]. Natural gas hydrates (NGHs) are a promising energy source with relatively low car-bon emissions and abundant global reserves that can provide solutions to reduce dependence on traditional fossil fuels and address energy shortages. However, the development of NGH faces challenges of technological complexity, high extraction costs, and environmental risks. The future economic potential depends on technological breakthroughs and the feasibility of large-scale exploitation, while the environmental sustainability of the extraction process must be ensured. In recent years, China has conducted several rounds of exploration and two successful experimental exploitations of buried hydrate areas in the South China Sea. According to the projections of the China Geological Survey, the South China Sea is rich in hydrate resources. As recorded in China Mineral Resources 2018, the gas hydrate resources in China’s waters are about 80 billion tons of oil equivalent, with the northern part of the South China Sea resource being 74.4 billion tons of oil equivalent [2,3]. How to develop marine gas hydrate efficiently and safely becomes one of the key points to ensuring energy security. The hydrate development methods that have been proposed can be broadly classified as (1) depressurization [4], (2) thermal stimulation [5], (3) CO₂ replacement [6], (4) chemical inhibitor injection, and (5) solid-state fluidization. Among these methods, depressurization and thermal stimulation are considered more feasible and have been used extensively in successful trial production tests worldwide. For example, as the more widely studied geological regions, the Nankai Trough in Japan [7,8,9] and the Shenhu Sea in China used depressurization and obtained stable gas production for a longer period [10,11,12]. However, a low-permeability gas hydrate reservoir hinders the propagation of a pressure drop during extraction and prevents efficient fluid transport from the formation to the wellbore, which in turn affects NGH decomposition and makes it difficult to enhance production [13,14].

As a key technology used to improve the reservoir permeability environment in low-permeability reservoirs, hydraulic fracturing has been largely researched and applied in the development of shale gas reservoirs as well as tight gas reservoirs [15,16,17]. The hydraulic fracturing method, which creates complex fractures through multi-stage fracturing segments to maximize the opening of new flow channels and obtain the highest oil recovery, has proven to be remarkably effective in its development for ultra-low-permeability reservoirs [18,19,20]. Hydraulic fracturing can form fractures in the reservoir to significantly increase the reservoir fluid flow area, which can be seen as a promising method for hydrate development when combined with conventional development methods such as depressurization. However, since hydrate development is still in the test state, there is no investigation of hydrate formation fracturing development at the engineering level, so whether hydraulic fracturing can effectively enhance the development production of NGH in the sea is still a key issue to be explored. The propagation of fractures in hydrate reservoirs and the effect of fractures on production enhancement have been explored by many scholars [21,22,23], and due to the difficulty of implementing large-scale exploitation experiments, the laboratory scale experiments as well as numerical simulations have become the choice of many researchers.

Since the decomposition of hydrates in the formation involves changes in several physical fields [24,25], most of the work on hydrate development simulation studies is performed by numerical simulators, among which TOUGH + HYDRATE, CMG-STARS, MH21-HYDRATES, and Hydrate-ResSim are extensively used by researchers for their superior performance. Sun et al. [26] studied the gas production potential of vertical fractures as well as horizontal fractures for the pressure-reducing extraction of horizontal wells using TOUGH + HYDRATE. Zhong et al. [11] investigated the production enhancement effect of a fracture network on hydrate reservoir development in the Shenhu area of site SH7 using an injection and extraction well pattern and focused on key fracture properties such as fracture network conductivity, horizontal or vertical branching fractures, etc. Xu et al. [27] used CMG-STARS to develop a simulation model for combined hydraulic fracturing in the depressurization extraction of three types of conventional hydrate reservoirs with different fracture parameters and investigated the gas production behavior, including fracture conductivity, location, size, and length. Liu et al. [28] proposed a method for the geothermal energy-assisted extraction of NGH, using Hydrate-ResSim to investigate factors such as the thermal recovery rate and hydrate-bearing layer (HBL) development performance, and to explore the technical feasibility of the scheme compared with conventional depressurization method extraction. In addition to the above numerical simulators, some scholars have used the COMSOL platform to perform multi-physical field coupling to simulate the decomposition behavior of hydrate in a reservoir, while Li et al. [29] proposed a novel modified method that combined hydraulic fracturing and depressurization assisted with sealing burdens, which analyzed the effects of factors such as radial length and permeability of the fracture zone on hydrate decomposition and gas production. However, since most of the above simulators do not have fracture description methods for fractures in hydrate reservoirs, currently, most of the production research on hydrate after fracturing is generally performed by local grid refinement (LGR) with increasing the permeability of that part of the refinement grid to simulate high-inflow channels [30,31]. Yu et al. [32] simulated horizontal fractures in the HBL by adding a thin horizontal layer of 0.01 m to the entire reservoir model, while Feng et al. [33] used an elliptical radially high-permeability zone to introduce fracture zones in the HBL. The generation of local grid refinement ensures accuracy for some simple fracture descriptions, but it is often difficult to simulate or the amount of computation increases significantly when the fractures show a trend of complexity. The method of fracture modeling has become a hot research topic in recent years [34,35,36]. Some classical models such as the equivalent permeability model, discrete fracture model (DFM), and embedded discrete model (EDFM) have been proposed [37,38,39] and many researchers have investigated their implementation [40,41] and adaptation refinement [42,43] in reservoir numerical simulators. However, these models have been rarely addressed in the research field of hydrate fracturing production research. Therefore, it is necessary to include the discussion of fracture models in the simulation studies of hydrate fracturing to increase production, which is useful to study the effect of increasing production in the actual fractures that present complications in the formation. For example, in a relatively new study, Liu et al. [44] developed a DFM-based HBL hydraulic fracturing-assisted pressure reduction exploitation model based on the fracture–matrix transfer relationship to achieve an accurate characterization of complex fracture shapes using a perpendicular bisection (PEBI) grid.

Among the hydrate production test cases around the world, China has made the most progress in exploitation, with two successful offshore test experiences to date [45,46]. In 2017 was the first trial recovery by the vertical well pressure reduction method for 60 days at the production test site in Shenhu Sea, South China Sea, in which the cumulative gas production reached 30.9 × 10⁴ m³ [47]. In 2020, the second offshore NGH production trial was carried out with success in the Shenhu area of the South China Sea at a depth of 1225 m [48]. This test production was carried out using the horizontal well pressure reduction method and achieved 30 days of stable gas production, with a cumulative gas production volume of 86.14 × 10⁴ m³ and an average daily gas production volume of 2.87 × 10⁴ m³ [49]. The breakthroughs of the two tests make China the first country in the world to successfully explore NGH in muddy silt sediments. However, its production level has not yet reached the standard of commercial exploitation, so it is necessary to investigate the effect of hydraulic fracturing to increase production in the Shenhu Sea. As shown in Figure 1, this paper aims to establish a new method to describe fractures in hydrate reservoir fracturing by implementing EDFM coupled with a conventional hydrate numerical simulator.

In this work, we firstly describe the process of constructing a hydrate production model and explain how we consider mechanisms such as multiphase flow; then, we validate the hydrate decomposition dynamics model used, and we also verify the feasibility by comparing the results of the EDFM model simulations with the LGR; and finally, based on all previous studies, we carry out an analysis of the influencing factors for three cases with different types of wells. A sensitivity analysis of key production factors such as fracture locations, fracture half-length, and fracture conductivity is performed in the results and discussion section by setting up three different combinations of hydraulic fracture-assisted depressurization for different well types. This study will provide useful references for future production enhancement and development of low-permeability hydrate reservoirs in marine areas.

2. Methodology

2.1. Geological Setting

The production test site in 2017 was located in the middle of the southeastern continental slope of the Shenhu area in the northern South China Sea between the Xisha Trough and the Dongsha Islands [47], and the hydrate reservoir in the sea around the test area SHSC-4 well was selected as the geological model for this study case.

From the analysis of logging data, the hydrate reservoir can be divided into three layers according to the distribution of different components in the pore space [50]: (1) The hydrate-bearing layer (HBL) located at 1495–1530 m has a distribution of hydrate and liquid water in the pore space. The porosity is 35% with a hydrate saturation of 0.34 and average permeability of 2.9 mD; (2) The three-phase layer (TPL), located at 1530–1545 m, has a pore space distribution of hydrate, liquid water, and free hydrocarbon gas. The porosity is 33%. The hydrate saturation is 0.31. The gas saturation is 0.078 [51]. The average permeability is 1.5 mD; (3) The free gas layer (FGL), located at 1545–1572 m, has pores filled with liquid water and free gas. The average permeability is 7.4 mD [47].

The size of the physical model is set to 280 × 280 × 137 m. As shown in Figure 2, the entire model was divided into five layers according to the components in the pore space, the overlying rock layer, HBL, TPL, FGL, and underlying rock layer, where the upper and lower layers are water-saturated, and the permeability and porosity are consistent with the parameters of their adjacent layers [50]. The simulation model is discretized as 56 × 56 × 48 in x, y, and z directions, where the x and y directions are uniformly discretized. The TPL is uniformly discretized into 15 meshes in z direction. The HBL is uniformly discretized into 18 meshes in z direction. The FGL is uniformly discretized into 9 meshes in z direction. The upper and lower rock layers are discretized into 3 meshes in z direction. The parameters in the numerical simulation model for hydrate development are shown in

Table 1. Some of the parameters used in the numerical simulations.

Parameters	Values
Rock density/(kg·m−3) [50]	2650
Rock thermal conductivity/(J·m−1·day−1·K−1) [52]	1.296 × 105
Thermal conductivity of hydrates/(J·m−1·day−1·K−1) [52]	3.395 × 104
Thermal conductivity of water/(J·m−1·day−1·K−1) [52]	5.183 × 104
Thermal conductivity of gas/(J·m−1·day−1·K−1) [52]	1.4 × 102
Geothermal gradients/(°C/km) [50]	48
Bottom of well pressure/(MPa) [53]	4.5
Sirw [53]	0.3
Sirg [53]	0.03

.

2.2. Gas Hydrate Equations

2.2.1. Kinetic Equations for Hydrate Decomposition

Some major hypotheses are proposed: a three-phase, three-component model is developed, where the three phases are water, gas, and hydrate phases (the hydrate is set as a solid phase), while the three components are divided into methane gas, water, and natural gas hydrate. The only gas produced by the decomposition of natural gas hydrates is methane, and the occurrence of decomposition reaction leads to changes in formation porosity, with the flow of liquid and gas phases in the pore space conforming to Darcy’s law. The numerical model includes mass conservation, energy conservation, and kinetic equations of hydrate decomposition.

In this study, the decomposition and formation of natural gas hydrate are considered as a reversible reaction of methane with water, expressed by the following equation:

$$C{H}_4(g)+n{H}_{2}O(l)\iff C{H}_4\cdot n{H}_{2}O(s)\pm heat,$$

(1)

where n is the hydration number associated with the hydrate structure and is taken as n = 5.75 [52,54]. The Kim–Bishnoi model [55] is generally used to describe the kinetics of hydrate decomposition by

$$m_{\mathrm{d}}={\lambda }_{\mathrm{d}}\left(\phi {\rho }_{\mathrm{w}}{S}_{\mathrm{w}}\right)\left(\phi {\rho }_{\mathrm{h}}{S}_{\mathrm{h}}\right)p_{\mathrm{e}}\exp \left(-{\displaystyle \frac{E}{RT}}\right)\left(1-{\displaystyle \frac{1}{K(p,T)}}\right),$$

(2)

where m_d is the intrinsic kinetic rates of NGH decomposition and formation; ${\lambda }_{\mathrm{d}}$ is the kinetic rate constants of hydrate decomposition and formation; $\phi$ is the intrinsic porosity, S is the phase saturation, and the subscripts g, w, and h represent gas, water, and hydrate, respectively; ${\rho }_{\mathrm{w}},{\rho }_{\mathrm{h}}$ are the densities of water and hydrate, respectively; p_e is the phase equilibrium pressure; E is the hydrate activation energy; R is the universal gas constant; p is the pressure; T is the temperature. K is the hydrate equilibrium or stability value, representing the effect of phase equilibrium during hydrate decomposition and formation [52], which can be expressed by the equation below,

$$K={\displaystyle \frac{k_1}{p}}\exp \left({\displaystyle \frac{k_2}{T-k_3}}\right),$$

(3)

where k₁, k₂, and k₃ are the fitted parameters.

2.2.2. Mass Balance Equation

The law of mass balance is used in the process of hydrate decomposition and formation reactions. The conservation of mass in this study can be expressed by the following equation [27].

$${\displaystyle \frac{\partial (\phi {\rho }_{\mathrm{g}}{S}_{\mathrm{g}})}{\partial t}}=\nabla \cdot ({\rho }_{\mathrm{g}}{v}_{\mathrm{g}})+{\dot{m}}_{\mathrm{g}}+q_{\mathrm{g}},$$

(4)

$${\displaystyle \frac{\partial (\phi {\rho }_{\mathrm{w}}{S}_{\mathrm{w}})}{\partial t}}=\nabla \cdot ({\rho }_{\mathrm{w}}{v}_{\mathrm{w}})+{\dot{m}}_{\mathrm{w}}+q_{\mathrm{w}},$$

(5)

$${\displaystyle \frac{\partial (\phi {\rho }_{\mathrm{h}}{S}_{\mathrm{h}})}{\partial t}}={\dot{m}}_{\mathrm{h}},$$

(6)

where v_g and v_w are the flow rate; $\dot{m}$ denotes the mass change due to hydrate dissociation; and q_g and q_w are the gas production rate and water production rate of the well, respectively.

2.2.3. Energy Conservation Equation

The law of energy conservation is observed in the process of hydrate decomposition and formation reactions [27]. The conservation of energy can be expressed by the following equation.

$$\begin{array}{l}\nabla \cdot \left({\lambda }_{\mathrm{c}}\nabla T\right)-\nabla \cdot \left({\rho }_{\mathrm{g}}{v}_{\mathrm{gr}}{H}_{\mathrm{g}}+{\rho }_{\mathrm{w}}{v}_{\mathrm{wr}}{H}_{\mathrm{w}}\right)+q_{\mathrm{g}}{H}_{\mathrm{g}}+q_{\mathrm{w}}{H}_{\mathrm{w}}+Q_{\mathrm{h}}=\\ {\displaystyle \frac{\partial }{\partial t}}\left[(1-\phi ){\rho }_{\mathrm{s}}{H}_{\mathrm{s}}+\phi ({\rho }_{\mathrm{h}}{S}_{\mathrm{h}}{H}_{\mathrm{h}}+{\rho }_{\mathrm{w}}{S}_{\mathrm{w}}{H}_{\mathrm{w}}+{\rho }_{\mathrm{g}}{S}_{\mathrm{g}}{H}_{\mathrm{g}}\right]\end{array},$$

(7)

where H is the enthalpy of each phase; and Q_h is the heat of hydrate dissociation, and ${\lambda }_{\mathrm{c}}$ can be described by the volume averaging of the effective heat conduction:

$${\lambda }_{\mathrm{c}}={\lambda }_{\mathrm{s}}(1-\phi )+\phi ({\lambda }_{\mathrm{h}}{H}_{\mathrm{h}}+{\lambda }_{\mathrm{g}}{H}_{\mathrm{g}}+{\lambda }_{\mathrm{w}}{H}_{\mathrm{w}}),$$

(8)

where ${\lambda }_{\mathrm{h}}, {\lambda }_{\mathrm{w}}, {\lambda }_{\mathrm{g}}$ are the thermal conductivity of each phase.

2.2.4. Relationship Between Porosity and Permeability

The effective permeability changes with porosity as gas hydrate decomposition occurs in the pore space [27,29]. The Carmen–Kozeny model is used to describe the relationship between effective permeability and porosity and can be expressed as

$$k=k_0{\left({\displaystyle \frac{\phi }{{\phi }_0}}\right)}^{\epsilon }{\left({\displaystyle \frac{1-{\phi }_0}{1-\phi }}\right)}^2,$$

(9)

where k is the effective permeability when the porosity is $\phi$, mD; and $\epsilon$ is an empirical parameter, which is taken as 5 in this study. This model describes the permeability change, which can be transformed into a form similar to the Tokyo model [56] under certain conditions. To verify the accuracy of using the Carmen–Kozeny model in simulation, we defined the initial porosity as 0.4, the $\epsilon$ as 2, and the permeability correction factor of the Tokyo model as 3. We compared the hydrate saturation from 0 to 1 for the two models with normalized permeability, as shown in Figure 3. The results show that both models fit well under the given conditions, indicating that it is possible to describe the variation between hydrate permeability and porosity using the Carmen–Kozeny model.

Since the hydrate phase is considered an immobile phase, its relative permeability is considered to be 0. The relative permeability model of the aqueous phase and the gas phase refers to the results of Moridis et al. [57], and the Van Genuchten model [58] to calculate the capillary pressure; the specific equation can be expressed as:

$$K_{rw}={\left({\displaystyle \frac{S_w-S_{irw}}{1-S_{irw}}}\right)}^{n_w},$$

(10)

$$K_{rg}={\left({\displaystyle \frac{S_g-S_{irg}}{1-S_{irw}}}\right)}^{n_g},$$

(11)

$$P_c=-p_0{\left[{\left({\displaystyle \frac{S_{\mathrm{w}}^{\ast }-S_{irw}}{1-S_{irw}}}\right)}^{-{\displaystyle \frac{1}{m}}}-1\right]}^{(1-m)},S_{\mathrm{w}}^{\ast }={\displaystyle \frac{S_{\mathrm{w}}}{1-S_{\mathrm{h}}}},$$

(12)

where K_rw and K_rg are the relative permeabilities of water and gas phases, md, respectively; S_irw and S_irg are the saturation of irreducible water and irreducible gas, respectively; n_w and n_g are the permeability reduction exponents of water and gas phases, respectively; P_c is the capillary force between gas and water phases; p₀ is the initial value of capillary pressure; m is an empirical parameter; and m = 0.45 is taken in this study.

2.3. Embedded Discrete Fracture Model

The embedded discrete fracture model (EDFM) is used in this work to describe hydraulic fractures in hydrate reservoirs, and EDFM is an effective method to simulate fractures with complex geometries. By “embedding” the fractures as well-like sources in the matrix, the model avoids local grid encryption and uses orthogonal structured meshing to increase the computational speed. The coupling between EDFM and the hydrate decomposition kinetic model is established by constructing a series of non-neighborhood connections (NNCs) and varying the grid conductivity coefficients of the fracture computational domain [17,59], as shown in Figure 4.

To maintain the pore volume of the fracture segment, an effective pore volume needs to be assigned to the fracture grid [41], which can be expressed as

$${\phi }_f={\displaystyle \frac{S_{seg}{w}_f}{V_{\mathrm{b}}}},$$

(13)

where ${\phi }_f$ is the effective porosity of the fracture cell; S_seg is the fracture segment area perpendicular to the fracture aperture; w_f is the fracture aperture; and V_b is the bulk volume of the cell assigned for the fracture segment.

The types of non-neighborhood connections can be divided into three categories, and the connection of a fractured segment to the matrix grid it penetrates is the first category, and the conductivity coefficient for this type is defined as shown in the following equation:

$$T_{f-m}={\displaystyle \frac{A_{f-m}{k}_{f-m}}{d_{f-m}}},$$

(14)

$$d_{f-m}={\displaystyle \frac{{\displaystyle {\int }_V{x}_{n}dV}}{V}},$$

(15)

where T_f−m is the conductivity coefficient between the matrix and the fracture; A_f−m is the area of the fracture segment embedded in the matrix; d_f−m is the average normal distance from the matrix to the fracture segment; and k_f−m is the summed average of the matrix permeability and fracture permeability. V is the volume of the matrix grid, and x_n is the distance from the volumetric microelements in the matrix to the fracture plane.

The conductivity of the fracture grid between single fracture segments is considered the second type of NNC, with the conductivity coefficients described by a two-point flux approximation, similar to the approximate simplification proposed by Karimi-Fard et al. [60], which can be expressed as

$$T_{seg}={\displaystyle \frac{T_1{T}_2}{T_1+T_2}},$$

(16)

$$T_1={\displaystyle \frac{k_f{A}_c}{d_{seg1}}},T_2={\displaystyle \frac{k_f{A}_c}{d_{seg2}}},$$

(17)

where T_seg is the conductivity coefficient of the fracture grid, k_f is the fracture permeability, A_c is the area of the common surface of two adjacent fracture segments, and d_seg1 and d_seg2 are the distances from the center of fracture segment 1 and fracture segment 2 to the common surface, respectively.

The connection between intersecting fracture segments is considered as the third type NNC, and the mass exchange for fracture intersection is approximated by setting a conductivity coefficient at the intersection of fracture segments according to Monifar et al. [59], as shown in the following equation:

$$T_{\mathrm{int}}={\displaystyle \frac{T_1{T}_2}{T_1+T_2}},$$

(18)

$$T_1={\displaystyle \frac{k_{f1}{w}_{f1}{L}_{\mathrm{int}}}{d_{f1}}},T_2={\displaystyle \frac{k_{f2}{w}_{f2}{L}_{\mathrm{int}}}{d_{f2}}},$$

(19)

where T_int is the conductivity coefficient at the intersection of the fracture segments, L_int is the length of the fracture intersection line, and d_f1 and d_f2 are the weighted averages of fracture sub-segment 1 and fracture sub-segment 2 to the intersection line, respectively, which can be expressed as

$$d_{f1}={\displaystyle \frac{{\displaystyle \underset{s_1}{\int }{x}_{n}d{S}_1+{\displaystyle \underset{s_3}{\int }{x}_{n}d{S}_3}}}{S_1+S_3}},$$

(20)

$$d_{f2}={\displaystyle \frac{{\displaystyle \underset{s_2}{\int }{x}_{n}d{S}_2+{\displaystyle \underset{s_4}{\int }{x}_{n}d{S}_4}}}{S_2+S_4}},$$

(21)

Considering the case of hydraulic fractures intersecting the wellbore, according to Monifar et al., this relationship is implemented in EDFM in the form of a specified well index (WI), as shown in the following equation:

$$WI={\displaystyle \frac{2\pi k_f{w}_f}{\ln (\frac{r_e}{r_w})}},$$

(22)

$$r_e=0.14\sqrt{L_s^2+H_s^2},$$

(23)

where L_s is the length of the fracture section, H_s is the height of the fracture section, r_e is the effective radius, and r_w is the radius of the wellbore.

2.4. Solution

The development of natural gas hydrate involves complex physical processes such as gas–liquid–solid phase change, temperature change, gas–water two-phase seepage, reservoir skeleton deformation, etc. In this study, we treat the hydrate phase as a solid phase component, which has been verified by the results of many researchers [52,61].

The implementation of EDFM is described as the following steps:

(1) Add to the original geologic model a number of fracture grids (computational domains) that do not exchange the medium with the matrix grids and that describe the fractures “embedded” in the matrix. The blocking of the medium exchange can be achieved by changing the transmissibility between the grids or by placing the null grids between the matrix and the fracture grids.

(2) To simulate the actual fracture volume in the discrete fracture network, the porosity of the computational domain needs to be determined from the fracture grid and the ratio of the actual fracture segment volumes. Equation (13) is referred and the grid properties of the computational domain are modified.

(3) The feature distances are determined from the intersection between the actual fracture segments and the matrix grid, and a first type NNC list can be generated. The first type of NNC is realized by defining a non-neighborhood connection between the matrix grid and the computational domain.

(4) The second and third type NNCs are identified by the spatial relationships between the actual fracture segments in the matrix grid. If the computational domain grids are all connected, this can be achieved by a change in the grid transmissibility in a specific direction. If they are not connected, this can be realized with the NNC.

(5) Determine whether the actual fracture section intersects the wellbore and if so, define a specific well grid in the computational domain and obtain the actual single well production, where the properties of the well grid are set concerning Equations (22) and (23).

3. Model Validation

3.1. Hydrate Decomposition Kinetic Model Validation

In this work, the kinetic equations for the decomposition and generation of natural gas hydrate are coupled with the heat and mass transfer equations; to verify the validity of the model, the simulation results of the model are compared with the resultant curves of Li et al. [62]. The physical dimensions of the model are 240 m × 240 m × 10 m. One-quarter of the model is selected for the analysis according to its symmetry, with a grid uniformly divided into 15 × 15 × 10 in the x, y, and z directions. The model simulates the effect of pressure reduction for Class II hydrate reservoirs with HBL in the upper 8 m of the model and the bottom water layer in the lower 2 m. A constant bottom well pressure of 4 MPa is used. Figure 5 shows the results of comparing the simulation results of the model in this study with the gas production curve of Case 1 by Li et al. As shown in Figure 5, the simulation results of the two models fit with high accuracy, which can illustrate the feasibility of the model constructed in this study.

3.2. Feasibility Verification of Embedded Discrete Fracture Model

The embedded discrete fracture model was used to describe the hydraulic fractured fractures in the reservoir. To verify the feasibility of the EDFM coupled hydrate extraction model, fractures are formed in a local grid refinement (LGR) with the same parameter settings to compare the gas production curves and the evolution of the physical parameters of the EDFM case in this study [63,64,65].

In this subsection, for the validation of the newly developed simulation method, three classical test cases are presented, as shown in Figure 6. To improve the convergence and computational speed of the LGR model in the validation case, we appropriately broaden the width of the fracture grid while keeping the fracture conductivity F_c at a constant value. To verify the effectiveness of this setting, Case 1 is chosen as the study object. The F_c is defined as 75 mD∙m, based on which we explore the comparison of gas production results for 1000 days of simulation for fracture widths of 0.1 m (k_f = 750 mD), 0.01 m (k_f = 7500 mD), and 0.001 m (k_f = 75,000 mD), respectively. As shown in Figure 7a, the simulated results for all three fracture widths are in high accordance, which indicates that the operation of appropriately widening the fracture grid while keeping the F_c at a constant value is feasible in this part of the validation case.

Case 1: The reservoir dimensions are 560 m × 560 m × 137 m, uniformly discretized into 56 × 56 grids along the X and Y directions, three grids for the overlying and underlying rock layers in the Z direction, five grids for the HBL, five grids for the TPL, and three grids for the FGL, and the rest of the conditions are consistent with the settings in the first test production model of Shenhu area built in this study. The F_c is set to 75 mD∙m. The grid where the fracture intersects the wellbore is appropriately refined.

Case 2: The reservoir size is 560 m × 560 m × 137 m. The rest of the conditions in the geological model are set in the same way as in Case 1. In this case, a horizontal well is used for production, and five fractures are set in the model; the fractures are bi-wing fractures and pass through the whole TPL. As a reference scheme, the F_c is set to 200 mD∙m, and the grid where the fractures intersect the wellbore is appropriately refined.

Case 3: The reservoir size is 560 m × 560 m × 137 m. The remaining conditions in the geological model are the same as in Case 1. In this case, one injection and one production vertical well are used, and two vertical wells are set on both sides of the model diagonal. Two fractures are set in the model, and the fractures intersect each other in a “+” shape, passing through the whole TPL in the longitudinal direction. As a reference scheme, the F_c is set to 500 mD∙m.

Figure 7 shows the comparison curves of the two models, with the simulation results presented in the form of gas production rates and cumulative gas production. The fitted results of EDFM simulation and LGR simulation are highly compatible with different extraction methods (vertical well, horizontal well, and one injection and one production with two vertical wells) and different fracture arrangements (vertical wells with double-wing fractures, horizontal wells with multiple sets of fractures, and “+” intersecting fractures). The errors of the three cases were statistically analyzed, in which the mean relative error (MRE) of each component was calculated as shown in

Table 2. Statistical analysis of errors for case validation.

Validation Cases	Curve Category	MRE
Case 1	Cumulative gas production	1.44%
	Gas production rate	1.41%
Case 2	Cumulative gas production	1.14%
	Gas production rate	1.05%
Case 3	Cumulative gas production	0.13%
	Gas production rate	0.22%

.

To reflect the results of the visual comparison and to increase the persuasiveness of the model, Figure 8, Figure 9 and Figure 10 are placed as a comparative result of the physical field evolution (both selected as examples for the twelfth layer where the horizontal wellbore is located), from which the feasibility of the model in this study can also be visualized. This shows that the coupling of EDFM in the hydrate reservoir development model is successful and applicable to subsequent studies. For multi-stage fractured horizontal wells, where fracture complexity is high, EDFM offers an advantage due to its ability to model the interaction of multiple fractures along the wellbore without the need for fine grid resolution. This makes it ideal for simulating large-scale fracture networks and their contribution to overall well productivity. In practical applications, particularly in complex fractured reservoirs, a hybrid approach that combines EDFM for broader fracture networks and LGR for single fracture areas could provide the optimal balance between efficiency and accuracy.

4. Case Study

In this section, EDFM is used to construct a series of production models for different well types with fractures, and the effect of different parameters of the fracture on the production increase in these production models is analyzed. We consider that some reservoir stimulation regions can be formed around the fracture during hydraulic fracturing. Therefore, we set a relatively high-permeability stimulation reservoir volume (SRV) around the fracture by changing the permeability, and the permeability of these areas was set to 30 mD, as shown in Figure 11. As the permeability of SRVs in the field can be affected by a variety of factors, this permeability value is derived from empirical values. Depending on the type of well, the model can be classified into the production of vertical fractured wells, multi-stage fractured horizontal wells, and fractured spiral multilateral wells, where the spiral multilateral well is constructed by referring to the results of Mao et al. [66].

4.1. Vertical Well Fracturing Production

4.1.1. Effect of Fracture Arrangement Layer

In this subsection, seven sets of simulations were performed as a comparison by setting the fractures in different layers to investigate the most suitable location for fracture arrangement in vertical fracture exploitation. The fracture was set to run through the layer in which it is located in the longitudinal direction, as shown in

Table 3. The setting of the fracture locations in the simulation cases.

Simulation Cases	Fracture Arrangement Layer
Run 1	HBL
Run 2	FGL
Run 3	TPL
Run 4	HBL + TPL
Run 5	HBL + FGL
Run 6	TPL + FGL
Run 7	HBL + TPL + FGL

for the specific simulation category settings. It can be divided into three main categories; that is, the fractures are located in a single layer, double layers, or the entire hydrate-bearing sediments. The perforation points of the vertical well in the simulation cases were kept consistent, and the fracture arrangement is shown schematically in Figure 11. In these simulation cases, the L_f was set to 52.2 m, and the F_c was set to 100 mD∙m.

Figure 12 shows the gas production as well as water production curves for 3000 days of simulated runs under different fracture locations. Figure 12a depicts the variation in cumulative gas production (V_g) of the seven simulation groups. It can be seen that in the seven simulation cases, the final V_g of Run 2 (FGL) is 2.468 × 10⁷ m³, while the cumulative gas production of Run 4 (HBL + TPL) is the highest, at 4.305 × 10⁷ m³, and the V_g of Run 7 is close to Run 4, at 4.286 × 10⁷ m³. The V_g of the entire simulation group shows an increasing trend with the increase in the fracture longitudinal height (L_H), but it is not linearly correlated with L_H. This indicates that the location of the fractured well is also one of the factors affecting production. In the three cases of single-layer fracture exploitation, HBL had the best fracture production increase, while FGL had the lowest production, which was also reflected in the case of double-layer fracture exploitation, as evidenced by the fact that Run 5 consistently produced more gas than Run 4. In addition, based on the comparison of the final production of Run 1 (L_H = 36 m) and Run 6 (L_H = 42 m), it can be seen that, where the final V_g of TPL + FGL is 3.497 × 10⁷ m³ and that of HBL is 4.146 × 10⁷ m³, the final cumulative gas production of Run 6, which has a larger longitudinal height of the fracture, is yet lower than that of Run 1, which can also indicate the greater production potential when HBL is chosen as the location for the fractured well.

Figure 12b shows the variation curve of the gas production rate (Q_g) in the simulated 3000-day production. For single-layer fracturing (Run 1–Run 3), the trend of Q_g can be divided into two stages. In the first stage, Q_g increases and reaches a peak, which occurs due to the existence of the drawdown between the reservoir and the wellbore, and due to the presence of a fracture that facilitates the propagation of the pressure drop, which provides a highly efficient pathway for the transportation of gas. In the second stage, the Q_g reaches a peak and then gradually decreases, which occurs because the decomposition rate of hydrate decreases due to the output of reservoir fluids that reduce the drawdown pressure. In contrast, for double-layer as well as three-layer fracturing, Q_g rises sharply in the pre-production period and then declines to some extent, which is particularly evident when HBL is involved (Run 4, Run 5, and Run 7), and it can be hypothesized that fracturing HBL can result in a higher Q_g in the production period.

Figure 12c, as well as Figure 12d, depicts the water production behavior of fracturing in different layers. It can be seen that in single-layer fracturing production, Run 3 has the lowest cumulative water production (V_w) of 3.823 × 10⁵ m³; Run 2 has the highest V_w of 4.007 × 10⁵ m³. Unlike the production dominance of HBL in V_g, the V_w of HBL single-layer fracturing is 3.832 × 10⁵ m³, which is extremely close to the simulation results of Run 3, which indicates that Run 1 has the highest gas-to-water ratio and the best production efficiency. In double-layer fracture exploitation, it can be found that Run 6 has the highest V_w and exceeds Run 7. It can also be seen in the variation in the rate of water production (Q_w) that Run 6 in the double-layer development remains high in the late stage of production, which is a disadvantageous factor for long-term production. Therefore, in combination with the variation in V_g, it can be judged that Run 6 has the lowest development economics.

Among seven simulation cases, Run 4 (HBL + TPL) and Run 7 achieved the highest V_g, while single-layer HBL fracturing also showed the best gas production, and FGL the worst. In multi-layer fracturing, HBL involvement resulted in higher gas production rates (Q_g) during early production, but Q_g declined as reservoir fluids reduced the pressure drawdown. Water production (V_w) showed that single-layer HBL fracturing had a high gas-to-water ratio and good production efficiency. However, in double-layer fractures, Run 6 had the highest water production and poor long-term economics due to sustained high water production rates. These results indicate that the fracture location, particularly targeting HBL, significantly affects gas recovery, with higher production potential and efficiency in certain fracture configurations.

4.1.2. Effect of Fracture Conductivity

This subsection discusses the effect of fracture conductivity (F_c) on the production increase in vertical well fracturing in hydrate reservoirs in the South China Sea test area. HBL was selected as the layer to be fractured, the L_f was set to 52.5 m, and the F_c was set to 10, 50, 100, 200, and 300 mD∙m, respectively.

Figure 13a, as well as Figure 13b, shows the gas production behavior for different fracture conductivity conditions. It can be seen that the presence of the fracture greatly increases the gas production. In the case of the change in V_g due to an increase in the same scale of conductivity, we can find that the V_g simulation results are 2.650 × 10⁷ m³ for F_c = 10 mD∙m and 3.923 × 10⁷ m³ for F_c = 50 mD∙m. Furthermore, the final production of 4.143 × 10⁷ m³ for F_c = 100 mD∙m provides only a small margin increase in V_g for a one-fold increase in F_c compared to the production of F_c = 50 mD∙m. From the above characteristics, it is not difficult to conclude that ever higher F_c means higher V_g, but the production does not increase linearly with the conductivity. It can be anticipated that after a certain value is reached, the increase in F_c will no longer be beneficial for the development of hydrate reservoirs in the South China Sea test area, which is also consistent with the findings of Zhong et al. [31].

From the variation in the gas production rate (Q_g) in Figure 13b, the overall curve variation can still be divided into two stages; specifically, the time of occurrence of peak Q_g is accelerated with the increment of F_c, which is due to the high-conductivity fractures providing more efficient fluid transport channels as well as the pressure drop propagation rate, which makes the hydrate decomposition rate faster. With the further growth of F_c, we found that a new stage of Q_g gradually emerged, a rapid increase in Q_g followed by a small decrease in the early stage of production for the 200 and 300 mD∙m simulation cases, which was caused by the high conductivity fractures accelerating the decomposition of hydrates around the fractures. In addition, the Q_g decreases significantly with the enhancement in F_c in the late stage of the simulation, where the final Q_g in the 300 mD∙m case is only 7032.1 m³/day, indicating that the effect of the high conductivity of the fracture to increase the production in the late production stage is not obvious, and can provide a better production increase only in the early stage. This is similar to the simulation results of Feng et al. [33].

Figure 13c as well as Figure 13d show the curves of water production behavior for different F_c conditions with V_w simulated results of 3.777 × 10⁵ m³, 3.824 × 10⁵ m³, 3.842 × 10⁵ m³, and 3.845 × 10⁵ m³ for F_c = 50, 100, 200, and 300 mD∙m, respectively. The results show that the Q_w increases with the increase in F_c, and the rate of growth of Q_w is fast in the early stage of development, but the increment rate decreases with time. In late production, the production pressure difference is reduced due to the continuous extraction of formation fluids, and the presence of high conductivity fractures exacerbates this process, resulting in a certain degree of decline in both V_w and Q_w in the later period.

4.1.3. Effect of Fracture Half-Length

This subsection discusses the effect of fracture half-length on the production increase in vertical well fracturing in hydrate reservoirs in the South China Sea test area. HBL was selected as the layer to be fractured, and the F_c was set to 100 mD∙m, while the fracture half-lengths were set to L_f = 0, 27.5, 52.5, 77.5, and 102.5 m, respectively.

Figure 14a, as well as Figure 14b, depicts the gas production behavior curves for 3000 days of simulation under different L_f conditions. From the variation in V_g in Figure 14a, compared with the simulation results (1.133 × 10⁷ m³) in the group without fractures (L_f =0 m), the V_g in the simulation cases with fractures all show a significant increase, with the final V_g for L_f = 102.5 m reaching 4.440 × 10⁷ m³, an increase of about 3.92 times. The reason is that the extension of the L_f expands the contact area between the reservoir and the fracture, which increases the pressure drop propagation and widens the area of hydrate decomposition. This is corroborated by the change in Q_g in Figure 14b. The time to peak Q_g shortens with increasing L_f, accelerating from 1220 days (L_f = 27.5 m) to 445 days (L_f = 102.5 m), and the peak also gains elevation with growing L_f, rising from 15,777.4 m³/day (L_f = 27.5 m) to 20,296.5 m³/day (L_f = 102.5 m). In addition, we observed a tendency for Q_g to decline with an increment in L_f in late production, which is like the effect of F_c. Combining the effects of fracture half-length and conductivity in the vertical well depressurization method reveals that the increase in production from the presence of fractures occurs mainly in the pre-production period.

Figure 14c,d depict the water production behavior curves for simulated 3000-day production under different L_f conditions. The enhancement in V_w, as well as Q_w, is not significant with increasing L_f, where V_w is very close for L_f = 52.5, 77.5, and 102.5 m, which are 3.824 × 10⁵ m³, 3.791 × 10⁵ m³, and 3.753 × 10⁵ m³. The final simulation results of Q_w are 18.29 m³/day, 13.20 m³/day, and 10.18 m³/day, which indicate that a higher L_f can reduce the water production and improve the gas–water ratio for long-term production and that increasing the L_f can improve the economy of hydrate reservoir development in the South China Sea to a certain extent.

4.2. Multi-Stage Fractured Horizontal Well Production

4.2.1. Impact of Multi-Stage Fracturing Horizontal Well Placement Layer

This subsection explores the effect of the arrangement of multi-stage fractured horizontal wells in each formation on production. The perforation points and fracture settings of the horizontal wells are shown in Figure 11, where the number of fractures is 7, the fracture spacing is 25 m, the L_f is 27.5 m, and the F_c is 100 mD∙m. The horizontal well sections are arranged in HBL, TPL, and FGL, respectively, and the fractures are set longitudinally through the whole layer where the horizontal wells are located.

Figure 15a,b depict the gas production behavior curves of horizontal well sections arranged in different layers. From the variation in V_g with time in Figure 15a, it is obvious that the V_g obtained by arranging the multi-stage fractured horizontal wells in HBL is the maximum, 5.444 × 10⁷ m³, and the V_g in FGL is the minimum, 3.802 × 10⁷ m³. The reason for this situation is that the perforation points and hydraulic fractures are all located in the horizontal section, whereas owing to the low permeability of the reservoir, the multi-stage fractured horizontal wells in the FGL cannot spread the pressure drop smoothly to the HBL and TPL during the process of pressure reduction, which makes the hydrate decomposition rate in these two formations low, and the decomposed fluid cannot be transported to the wellbore in time. Overall, multi-stage fractured horizontal wells have a greater yield advantage (final V_g of 5.444 × 10⁷ m³ for HBL) compared to vertical wells (V_g of up to 4.305 × 10⁷ m³ for Run 4).

Figure 15b shows the curve of Q_g versus time. HBL as well as TPL still conform to the two stages of gas production rate variation, and the former reaches its peak faster than the latter and has a higher peak than the latter. However, when the horizontal well section is arranged in the FGL, Q_g displays three phases: the first stage (0~56 days), in which Q_g first decreases and then rises, reaches a peak, and then falls to a stable gas production stage; the second stage (56~620 days), which is a stable gas production stage, and the gas production rate remains stable and gradually rises at about 530 days; and the third stage (620 to 3000 days), in which the curve shows a decreasing trend, corresponds to the second stage in the other two cases. The reason for the occurrence of these three stages is that the FGL has the highest matrix permeability of the three layers and contains free methane gas, and the gas flows into the wellbore more quickly depending on the fracture at the maximum production pressure difference in the pre-production period, so the peak gas production rate occurs earlier. After a period of production, the pressure drop gradually propagates to TPL and HBL owing to the fluid recovery in FGL and the influence of fractures, and the gas produced by hydrate decomposition is transported to the horizontal section of the wellbore in FGL, leading to the stage of stable gas production, after which the rate of hydrate decomposition slows down due to the decrease in the production pressure difference, so the third stage, i.e., the process of a gradual decrease in the gas production rate, occurs.

Figure 15c, as well as Figure 15d, depicts the water production behavior curves of horizontal well sections arranged in different layers, where HBL < TPL < FGL in the comparison of V_w with 2.971 × 10⁵ m³, 3.504 × 10⁵ m³, and 3.861 × 10⁵ m³, respectively. Q_w presents a trend of FGL < HBL < TPL in the early period (0~530 days), whereas after 530 days, the Q_w of HBL will be the lowest and has been showing a gradual decrease, and the final simulation result is only 5.032 m³/day, which has a further decrease compared with the straight wells (Run 1), which indicates that the exploitation economy of multi-stage fractured horizontal wells is better than that of vertical wells with single-layer fracture development.

Figure 16, Figure 17 and Figure 18 show the images of physical field evolution with time for different fracture arrangement layers in multi-stage fractured horizontal wells, respectively. In Figure 18, the image of hydrate concentration evolution, there are additional areas of undecomposed hydrate present in the simulated case of FGL fracturing. Combining the evolution of temperature and pressure, it can be found that the propagation of pressure waves is more difficult, and the overall low-pressure region is mostly restricted to the FGL, for which the influence of the hydrate decomposition region exists only in the lower formation of the TPL. When the horizontal well section and the fracture are arranged in the HBL, the presence of the fracture connects the upper boundary of the TPL to the entire HBL, resulting in the largest area of low pressure and low temperature in the formation around the wellbore, which is particularly evident in the 400-day comparison. In addition, in the observation of hydrate concentration, it can be found that hydrate decomposition is the most adequate at the upper and lower boundaries of HBL, and most of the hydrate in TPL is decomposed, which also indicates that arranging the fracture and horizontal well section in HBL is the optimal choice.

In the evolution of the temperature field in Figure 17, the low-pressure and low-temperature regions are significantly less in the FGL simulation case than in the other two cases during the same time. There is a more obvious replenishment of high-temperature fluids from the bottom layer to the upper layer in the TPL case, but the effect on the upper HBL of the wellbore is smaller, so there are still a large number of incompletely decomposed regions in the HBL in the final hydrate concentration comparison.

4.2.2. Effect of Fracture Conductivity in the Multi-Stage Fractured Horizontal Well

In this subsection, the effect of different F_c of fractures on production is investigated in a multi-stage fractured horizontal well production, where the horizontal well section is arranged in HBL, the fractures are set longitudinally through HBL, the number of fractures is 7, the fracture spacing is 25 m, the fracture half-length is set to L_f = 27.5 m, and the conductivity is F_c = 10, 50, 100, 200, and 300 mD∙m.

Figure 19a, as well as Figure 19b, depicts the gas production behavior curves for different F_c in a multi-stage fractured horizontal well production. Compared to the V_g curves for different F_c in a vertical fractured well production in Figure 13a, the gas production obtained for the same value of F_c in a multi-stage fractured horizontal well is higher, where the V_g for Fc = 300 mD∙m is 5.522 × 10⁷ m³, which is approximately a 27% increase compared to a vertical fractured well (V_g of 4.341 × 10⁷ m³ for vertical fracturing). In addition, the flattening of the production curve is more pronounced with increasing F_c in multi-stage fractured horizontal wells, and the curves in the figure have smaller gaps at later stages, indicating that fracture conductivity is an important factor affecting horizontal well production, but the benefits from significantly increasing this indicator are smaller.

From the variation in Q_g in Figure 19b, it is simple to find that the higher F_c in the production of multi-stage fractured horizontal wells has a more significant effect on the reduction in Q_g caused by the later stages of production compared to vertical fracturing, and this phenomenon appears earlier. Compared to the rate curve reversal in Figure 13b, which occurs around 1940 days, this happens sooner in multi-stage fractured horizontal wells to around 1100 days, and the Q_g of F_c = 200 and 300 mD∙m is even lower than the simulated group without fractures at the end of the simulated 3000 days, which indicates that in the later stages of the development of multi-stage fractured horizontal wells, the increase in fracture conductivity reflects more of a negative effect on production.

Figure 19c, as well as Figure 19d, depicts the water production behavior curves of the multi-stage fractured horizontal wells under different F_c conditions; compared to the gas production behavior, the changes caused by the elevated F_c in the V_w as well as Q_w curves remain at a more negative level, where the water production has tended to be consistent at 3000 days for the four simulated cases where the F_c increases from 50 to 300 mD∙m, with 2.925 × 10⁵ m³, 2.951 × 10⁵ m³, 2.970 × 10⁵ m³, and 2.979 × 10⁵ m³, respectively. In addition, for the F_c research, it can be found that the 10 mD∙m and no-fracture simulation groups are closer in the production dynamics of gas and water production, which suggests that in the case of horizontal wells, more perforation points provide more spheres of influence with the reservoir, which replaces the function of the low-conductivity fracture to a certain extent.

4.2.3. Effect of Fracture Half-Length in Multi-Stage Fractured Horizontal Wells

In this subsection, the effect of different fracture half-lengths on the production in a multi-stage fractured horizontal well is investigated; the rest of the conditions remain the same, changing the fracture half-length, which is set to L_f = 0, 12.5, 27.5, 37.5, and 52.5 m, respectively.

Figure 20a,b depict the gas production behavior curves of multi-stage fractured horizontal wells under different fracture half-length conditions. The change in L_f has a similar effect on changing the F_c for multi-stage fractured horizontal well production, with V_g of 3.674 × 10⁷ m³ for no fracture (L_f = 0 m). In contrast, the smallest fracture half-length (L_f = 12.5 m) in this study yielded a V_g of 5.369 × 10⁷ m³, which is an increase of about 46%. Moreover, it can also be judged from the closeness of the curves that the increase in L_f does not have a decisive influence on the change in production, as the characteristics of the change in the gas production rate also indicate. The decrease in the rate in the late production period as well as it being very pronounced in Figure 20b are detrimental to long-term gas production.

Figure 20c as well as Figure 20d depict the water production behavior curves of the multi-stage fractured horizontal wells under different L_f conditions. It can be observed that in Figure 20c, there is a phenomenon of lower V_w obtained for higher L_f cases, with a reversal of the magnitude of water production around 1000 days. In the final simulation results, V_w is 3.018 × 10⁵ m³ (L_f = 12.5 m), 2.971 × 10⁵ m³ (L_f = 27.5 m), 2.936 × 10⁵ m³ (L_f = 37.5 m), and 2.884 × 10⁵ m³ (L_f = 52.5 m), respectively. Combined with the Q_w curve, the increase in fracture half-length increases V_w slightly in the early stage of production but can suppress water production in the later stage.

4.3. Spiral Multilateral Well Fracturing Production

4.3.1. Effect of the Number of Fractured Branches

This section investigates the effect of the number of wellbore branches (N_f) on production in fractured wells in spiral multilateral wells (N_f is referred to as “fractured branches” in this section to distinguish it from normal wellbore branches).

In this section, the fractures in each branch are set to L_H = 6 m to reduce the effect of fracture vertical height, and since HBL has been found to have higher gas production potential in the study of fracture layers in vertical and horizontal wells, only the fractures in the wellbore branches in HBL are used in this section. The fracture spacing is 10 m, the number of fractures is 4, the F_c is 100 mD∙m, and the L_f is 12.5 m.

Figure 21a,b depict the gas production behavior curves for different numbers of fractured branches. Due to the presence of a certain number of wellbores and the distribution of numerous perforation points in each layer of the spiral multilateral wells, the impact area of the pressure drop is increased, which to some extent replaces the function of the artificial fracture section. Combined with the variation in V_g in Figure 21a, there is no more significant production enhancement after fracturing the area where some of the wellbore branches are located. In contrast to the case without fracturing (V_g is 4.889 × 10⁷ m³), the V_g is 4.974 × 10⁷ m³, 5.166 × 10⁷ m³, and 5.266 × 10⁷ m³ for fracturing branches N_f of 1, 3, and 5, respectively, and the maximum increase is only about 7.7%.

Figure 21c, as well as Figure 21d, depicts the water production behavior curves for a different number of fracture branches, like the trend presented in Figure 20c in the variation in L_f in the studied horizontal wells, which shows a reversal of V_w around 1730 days, with the simulated water production results of no fractures overtaking the fracture branches to reach the maximum in the late production period. The V_w is 3.336 × 10⁵ m³ for no fracturing and 3.182 × 10⁵ m³ for N_f = 5, which shows that fractures are conducive to improving increased water production in the later stages of production.

4.3.2. Effect of Fracture Conductivity

Figure 22a,b depict the effect of different fracture conductivity on gas production behavior in fractured spiral multilateral wells, where the fracture branches are set to 3 and F_c is set to 10, 100, and 300 mD∙m, respectively, and the rest of the parameters are kept constant. At the end of the 3000-day simulation, V_g was 5.214 × 10⁷ m³ for F_c = 300 mD∙m, and 4.889 × 10⁷ m³ in the simulation for the no-fracture group, an enhancement of only about 6.6%. From the variation in V_g with time, it can be seen that changing the fracture conductivity has improved the production of spiral multilateral wells. Still, the improvement in production by F_c is much lower compared to vertical and horizontal wells. In particular, in both water and gas production curves in Figure 22, the curve with F_c = 10 mD∙m almost coincides with the curve of the simulation result without fracturing, which also occurs in the simulation of the horizontal wells but is more obvious in the spiral multilateral wells, and it can be seen that, for the more complicated well types, the low-conducting fracture is no longer effective in increasing the production.

4.3.3. Effect of Fracture Half-Length

Figure 23 depicts the effect of different fracture half-lengths on gas production and water production behavior in fractured spiral multilateral wells with F_c set at 100 mD∙m and all parameters set as constant except L_f. There is a certain amount of an increase in V_g in the simulation results for different L_f but raising L_f did not have a significant increase in production. At the end of the 3000-day simulation, the V_g of Lf from 12.5 m to 52.5 m was 5.166 × 10⁷ m³, 5.213 × 10⁷ m³, 5.238 × 10⁷ m³, and 5.280 × 10⁷ m³, respectively, with a maximum incremental production of V_g of only about 2.2%. In addition, in the study of the water production behavior, it can be found that all the simulations are extremely close to each other in the curves, but there is a certain degree of reduction relative to the simulation group without fracturing; however, this effect does not imply that there is a need to increase the L_f deliberately, and in the case of the synthesis, the L_f is an insensitive factor for the spiral multilateral wells, and it is not very helpful for the increase in the production.

4.4. Impact of Stimulated Reservoir Volume Size

With the horizontal well as the research object, this subsection analyzes the effect of different stimulated reservoir volume sizes on production enhancement. The horizontal well section is arranged in HBL, and the conductivity is F_c = 100 mD∙m. Four simulation groups are set up as SRV-S-1, SRV-S-2, SRV-S-3, and SRV-S-4, where the stimulated size is set as shown in Figure 24. SRV-S-1 is for a single fracture, i.e., no stimulated volume; SRV-S-2 indicates that the stimulated size is limited to the grid where the main fracture is located; SRV-S-3 indicates that the stimulated size includes the area of one grid outside the fractures; and SRV-S-4 indicates the area of two grids outside the grid of the main fracture.

Figure 25 shows the dynamics of gas as well as water production for the simulated 3000 days at different stimulated reservoir volume sizes, with V_g of 5.092 × 10⁷ m³, 5.404 × 10⁷ m³, 5.444 × 10⁷ m³, and 5.471 × 10⁷ m³, and an increase in production of about 6.13%, 0.74%, and 0.50%, respectively. It is evident that as the stimulated reservoir volume size increases, V_g increases more only between SRV-S-1 and SRV-S-2. The change in the curve of Q_g suggests that larger stimulated reservoir volume sizes will allow the peak to advance but the decline to be more pronounced at a later stage. V_w are 2.951 × 10⁵ m³, 2.987 × 10⁵ m³, 2.971 × 10⁵ m³, and 2.948 × 10⁵ m³, respectively. It is not difficult to see from the visual change in the curves that, similar to V_g, V_w also has a more obvious change only between SRV-S-1 and SRV-S-2, but the water production of the four cases is close to unanimous in 3000 days, which indicates that the water production will reach the peak value earlier when stimulated reservoir volume is considered. The variation in Q_w also indicates that the case considering the stimulated reservoir volume has a high rate of water production in the early stages, but declines significantly in the later stages of production.

Figure 26 shows schematic diagrams of the evolution of the pressure field, temperature field, and hydrate concentration field after 3000 days of simulation, and it can be seen from the variations that the overall pressure, temperature, and hydrate saturation of the reservoirs in the remaining three cases decreased compared to SRV-S-1, which had no stimulated reservoir volume. It is not difficult to see that when the stimulated reservoir volume is enlarged, the complete decomposition area of hydrate at the perforation point is obviously diminished, and the decomposition trend of hydrate within the stimulated reservoir volume also decreases, and after synthesizing the results of gas production curves, it can be deduced that too large a stimulated reservoir volume size does not effectively promote the decomposition of hydrate. Only when the V_g and Q_g between SRV-S-1 and SRV-S-2 changed greatly, the further expansion of the transformation scale, i.e., the change between SRV-S-2, SRV-S-3, and SRV-S-4, was not obvious, which suggests that the increase in the production of hydrate reservoirs does not require the size of stimulated reservoir volume, and that it is only necessary to form the stimulated reservoir volume size of about 2–8 m around the main fracture in order to obtain a higher production efficiency and increase in production.

In this section, we analyzed the gas production performance across the three types of hydraulically fractured wells; the most substantial production increase was observed in multi-stage fractured horizontal wells, while the spiral multilateral wells exhibited the smallest production change. Considering production rates, horizontal wells emerge as the optimal choice when hydraulic fracturing is employed. The presence of high-conductivity fractures accelerates hydrate decomposition around the fracture, advancing the peak gas and water production rates as the fracture conductivity increases. However, a notable decline in production is observed during the later stages of simulation, suggesting that the positive impact of fracture conductivity on production diminishes over time. This indicates that high-conductivity fractures primarily enhance production during the early phases of depressurization. Additionally, under the condition that there is no interference between SRVs, it is believed that a moderate increase in the stimulated reservoir volume (SRV) can effectively enhance hydrate decomposition.

5. Conclusions

In this paper, by coupling EDFM with the NGH formation to study the production increase effect of hydraulic fracturing-assisted depressurization development, we set up a case study of fracturing production in different well types after considering the stimulated reservoir volume, and analyze the gas and water production behaviors under different fracture locations, fracture conductivity, and fracture half-length conditions. Based on the results of this study, a few conclusions were drawn as follows:

(1) For the geological conditions in the Shenhu first test area, considering the combination of gas and water production, the HBL is the most economical layer to develop even if multiple layers are hydraulically fractured, and the HBL is also the optimal horizontal wellbore location for multi-stage fractured horizontal wells. In the case of vertical and horizontal fractured wells, fracture conductivity is an important factor in production, whereby greater F_c means higher gas production. However, the production does not increase linearly with the fracture conductivity.

(2) In the study of the gas production behavior of the three well types with hydraulic fracturing, the most significant increase in production was obtained in the multi-stage fractured horizontal wells, with the least change in production in the spiral multilateral wells. The optimal type choice is a horizontal well if hydraulic fracturing is used in production.

(3) The high-conductivity fracture accelerates the decomposition of hydrates around the fracture, and the time of peak gas and water production rates are accelerated with the increase in F_c, but the decline is also more and more obvious at the later stage of simulation. This indicates that the effect of the high conductivity of the fracture to improve the production is no longer apparent at the later stage of production, which can provide a better production increase only at the early stage.

(4) The presence of a certain number of wellbores and perforation points in each sub-stratum of the spiral multilateral wells facilitated the propagation of a pressure drop in the reservoir, replacing to some extent the effect of a fracture production increase, thus making the production gain from increasing F_c and L_f weakened.

(5) Appropriate stimulated reservoir volume size can obtain high gas production and production efficiency. In the study case, forming a stimulated reservoir size of about 2–8 m around the main fracture can achieve the effect of promoting the decomposition of hydrate.