This study considers a fractured vertical well in a two-dimension (2D) naturally fractured gas hydrate reservoir. In order to consider heat transfer and to characterize the complex fractures, we make the following assumptions to simulate gas-water flow in hydrate reservoirs.

(1) hydrate is treated as an immobile solid phase;

For the matrix system, the mass conservation equation of the gas and water can be expressed as (Hong and Pooladi-Darvish, 2003): ∇ ρ l ν l m + ρ l q l , m f + m l = ∂ ∂ t ρ l V v S l m , l = g , w (1)

For the intersecting fracture system, the mass conservation equation of the gas and water can be expressed as: ∇ ρ l v l f + ρ l q l , w + ρ l q l , m f + ρ l q l , f f + m l = ∂ ∂ t ρ l V v S l f , l = g , w (2) where ν _l is volumetric flow rates, which can be expressed as: v l = A c k k r l μ l ∇ P l (3)

For hydrate, the mass conservation equations of the matrix and fracture systems are (Hong and Pooladi-Darvish, 2003): − m h m , f = ∂ ∂ t ρ h V v S h m , f (4) In Eqs 1–4, subscript l represents gas or water; subscript h represents hydrate; the superscripts m and f represent the matrix and the fracture, respectively; k represents absolute permeability; k _rl represents fluid relative permeability; P represents pressure; μ _l represents viscosity; ρ _l represents density; m _l represents the mass accumulation of gas and water due to hydrate dissociation; m _h represents the mass of hydrate dissociation; V _v represent void pore volume; S _l represents gas and water saturation; S _h represents hydrate saturation; q _l,mf represents volumetric flow rates between matrix and fracture; q _l,ff represents volumetric flow rates between fracture and fracture; q _l,w represents volumetric flow rates between fracture and well.

For the EDFM method, unrefined and structured grids are generally used to characterize complex natural fractures, and the finite difference method can be used as the numerical method to discretize mass conservation equations (Yan et al., 2014). For the right-hand-side of Eqs 1, 2, we can have the following discrete form: ∂ ∂ t ρ l V v S l = V v ρ l n + 1 c s + c l S l n + 1 P n + 1 − P n + V v ρ l n + 1 S l n + 1 − S l n d t (5) where c _s represents rock compossibility coefficient; c _l represents fluid compossibility coefficient.

The discrete form of fluid flow between matrix grids in Eq. 1 can be written as: ∇ ρ l ν l m = ρ l n + 1 T m , i j k r l n + 1 μ l n + 1 P j m − P i m n + 1 (6)

Likewise, the discrete form of fluid flow between fracture grids in Eq. 2 can be written as: ∇ ρ l ν l f = ρ l n + 1 T f , i j k r l n + 1 μ l n + 1 P j f − P i f n + 1 (7)

The discrete form of mass transfer between matrix grid and fracture grid in Eqs 1, 2 can be written as: ρ l q l , m f = ρ l n + 1 T m f k r l n + 1 μ l n + 1 P i m − P j f n + 1 (8)

The discrete form of mass transfer between intersecting fracture grids in Eq. 2 can be written as: ρ l q l , f f = ρ l n + 1 T f f , i j k r l n + 1 μ l n + 1 P j f − P i f n + 1 (9)

The discrete form of mass transfer between the wellbore and fracture grid in Eq. 2 can be written as: ρ l q l , w = ρ l n + 1 k r l n + 1 μ l n + 1 2 π k i Δ h ln r e / r w P i f − P w f n + 1 (10)

In Eq. 10, r _w is the wellbore radius, and r _e is the equivalent wellbore radius which can be calculated with the Peaceman formulation (Peaceman, 1978), and the details for calculating m _l in Eqs 1, 2 can be found in Supplementary Appendix SA.

In Eqs 6–9, the terms T _m , T _f , T _mf , and T _ff represent the transmissibility factors between different connections of matrix and fracture grids. Figure 1 presents the schematics of different connections of matrix and fracture grids. Three different types of non-neighbor connections (NNCs) in addition to the connections between neighboring matrix grids can be found in Figure 1, including 1) type #1 NNC, the connections between the embedded fracture segments and matrix grids; 2) type #2 NNC: the connections between neighboring fracture segments of a single fracture; and 3) type #3 NNC: the connections between intersected fractures segments. For the fluid flow between matrix grids, taking x-direction as an example, the transmissibility factor between neighboring matrix grid i and j can be written as: T m , i j = 2 T m , i T m , j T m , i + T m , j (11) where T m = k m Δ y Δ z Δ x (12) In Eq. 12, the k _m represents matrix grid permeability; ∆x represents x-direction matrix grid size; ∆y represents y-direction matrix grid size; ∆z represents z-direction matrix grid size.

For type #1 NNC, the matrix-fracture transmissibility factor can be expressed as (Lee et al., 2001): T m f = 2 A f k m d f m (13)

In Eq. 13, the A _f is the area of the cross section between matrix grid and fracture segment, and d _fm is the average normal distance between the embedded fracture segment and matrix grid. For type #2 NNC, the transmissibility factor of fluid flow between neighboring fracture segments of a single fracture is calculated as (Lee et al., 2001): T f , i j = 2 T f , i T f , j T f , i + T f , j , T f , i = k f i w i Δ h i Δ l f i (14)

In Eq. 14, the k _fi is permeability of ith fracture segment; ∆l _fi , ∆h _i , and w _i represent the length, height, and width of ith fracture segment, respectively. For type #3 NNC, the transmissibility factor of fluid flow between segment i of fracture 1 and segment j of fracture 2 can be expressed as (Lee et al., 2001): T f f , i j = T i , 1 T j , 2 T i , 1 + T j , 2 , T i , 1 = k f i , 1 w i Δ h d s e g 1 , T j , 2 = k f j , 2 w i Δ h d s e g 2 (15)

In Eq. 15, d _seg is the distance from the fracture segment to the intersection.

The residual mass conservation equations in terms of the finite difference method can be expressed as (Song and Liang, 2009): M a t r i x : R l m , i n + 1 = V v ρ l n + 1 c s + c l S l n + 1 P n + 1 − P n + V v ρ l n + 1 S l n + 1 − S l n i m d t − ∑ j = 1 n c ρ l n + 1 T m , i j k r l n + 1 μ l n + 1 P j m − P i m n + 1 − ∑ j = 1 n f ρ l n + 1 T m f k r l n + 1 μ l n + 1 P j f − P i m n + 1 − m l , i n + 1 F r a c t u r e : R l f , i n + 1 = V v ρ l n + 1 c s + c l S l n + 1 P n + 1 − P n + V v ρ l n + 1 S l n + 1 − S l n i f d t − ∑ j = 1 n c ρ l n + 1 T f , i j k r l n + 1 μ l n + 1 P j f − P i f n + 1 − ρ l n + 1 T m f k r l n + 1 μ l n + 1 P j m − P i f n + 1 − ∑ j = 1 n i n ρ l n + 1 T f f , i j k r l n + 1 μ l n + 1 P j f − P i f n + 1 − m l , i n + 1 (16)

For the hydrate, we can have: M a t r i x : R h m , i n + 1 = V v ρ h n + 1 c s S h n + 1 P n + 1 − P n + V v ρ h n + 1 S h n + 1 − S h n i m d t + m h , i n + 1 F r a c t u r e : R h f , i n + 1 = V v ρ h n + 1 c s S h n + 1 P n + 1 − P n + V v ρ h n + 1 S h n + 1 − S h n i f d t + m h , i n + 1 (17)

For the matrix system, the energy conservation equation can be expressed as (Hong and Pooladi-Darvish, 2003): ∂ U m ∂ t = ρ l ν l H l m + ∇ λ e f f ∇ T m + U m f + U d i s s (18) U m = V v U l S l ρ l + U h S h ρ h + U r V r m (19)

In addition, the energy conservation equation for the fracture system is: ∂ U f ∂ t = ρ l ν l H l f + ∇ λ e f f ∇ T f + U m f + U f f + U f w + U d i s s (20) U f = V v U l S l ρ l + U h S h ρ h + U r V r f (21) In Eqs 18–21, U ^m and U ^f represent internal energy of matrix and fracture system; U _mf represents energy exchange between matrix and fracture system; U _ff represents energy exchange between different fractures; U _fw represents energy exchange between fracture and well; U _diss represents heat absorbed by hydrate decomposition; U _l represents internal energy of gas and water; U _h represents internal energy of hydrate; U _r represents internal energy of rock; H _l represents enthalpy of gas and water; T represents temperature; λ _eff represents effective thermal conductivity; V _r represents rock volume.

The effective thermal conductivity λ _eff and enthalpy of the water and gas H _l are given as follows (Hong and Pooladi-Darvish, 2003; Esmaeilzadeh et al., 2008): λ e f f = ϕ v S w λ w + S g λ g + S h λ h + 1 − ϕ v λ r (22) H l = ∫ T r e f T l C p l d T (23) In Eqs 22, 23, ϕ _v represents void porosity; S _w , S _g , and S _h represent water, gas, and hydrate saturation respectively; λ _w , λ _g , and λ _h represent water, gas, and hydrate thermal conductivity respectively; T _ref represents reference temperature; C _pl represents fluid heat capacity.

The internal energy U _l is accumulated energy which neglects the mechanical work of a fluid phase. For gas and water, the relationship between U _l and H _l is (Esmaeilzadeh et al., 2008), U l = H l − P l ρ l (24)

In Eq. 24, P _l is phase pressure and ρ _l is density. Energy exchange U _mf between jth matrix and the ith fracture can be expressed as (Shi, 2021): U m f = ρ l q m f H l + C I λ m f T i f − T j m (25) In Eq. 25, CI is the connectivity index between the matrix and fractures (CI = A _f /d _fm ); λ _mf is the average thermal conductivity between the matrix and the fractures.

Energy exchange U _ff between ith segments of fracture 1 and the jth segments of fracture 2 can be expressed as (Shi, 2021): U f f = ρ l q f f H l + T I f f λ f f T i f 1 − T j f 2 (26)

In Eq. 26, λ _ff is the average thermal conductivity between different fractures; TI _ff is the transmissivity factor between different fractures which can be expressed as (Shi, 2021): T I f f = 2 T I f 1 T I f 2 T I f 1 + T I f 2 , T I f 1 = w f i Δ h Δ l f i (27)

Energy exchange U _fw between well and the ith fracture can be expressed as: U f w = ρ l q f w H l (28)

In general, the discrete form for the energy conservation equation of matrix and fracture grids can be expressed as (Song and Liang, 2009): Matrix : R t m , i = U n + 1 − U n i m Δ t − ∑ j = 1 n ρ l n + 1 ν l n + 1 c p l T n + 1 − T r m − ∑ j = 1 n λ m , j T j n + 1 − T i n + 1 m − ∑ j = 1 n f ρ l n + 1 q m f n + 1 c p l T j f − T r n + 1 + C I λ ¯ m f T j f − T i m n + 1 − m h , i n + 1 H h Fracture : R t f , i = U n + 1 − U n i f Δ t − ∑ j = 1 n ρ l n + 1 ν l n + 1 c p l T n + 1 − T r f − ∑ j = 1 n λ m , j T j n + 1 − T i n + 1 f − ρ l n + 1 q m f n + 1 c p l T j m − T r n + 1 − C I λ ¯ m f T j m − T i f n + 1 − ∑ j = 1 n i n ρ l n + 1 q f f n + 1 c p l T j f − T r n + 1 + T I f f λ ¯ f f T j f − T i f n + 1 − ρ l n + 1 q f w n + 1 c p l T i f − T r n + 1 − m h , i n + 1 H h (29)

For Eqs 16, 17, 29, values of k _rl , μ _l , ρ _l, and H _l are taken from the phase upstream region and can be solved with a fully implicit coupling (see Supplementary Appendix SB).

In this case, we established a 155×155×30 m³ hydrate reservoir with a single fracture to conduct the validation. Figure 2 gives a top view of the global grid system used to simulate the production of the established model. As shown in Figure 2, the local grid refinement (LGR) was used to simulate the hydraulic fracture. The reservoir domain is discretized into 31 × 31 = 961 matrix grids in the EDFM model. The bottom-hole pressure (BHP) of the fractured vertical gas well which locates at the center of the hydrate reservoir is set to be a constant value of 5 MPa. The fracture length is 35 m, and the fracture conductivity, which is 1,000 mD∙m, used in CMG STARS is consistent to that used in the EDFM model.

Figure 3 compares the calculated results from the proposed EDFM model to those from CMG STARS. Figures 3A,B illustrate the gas production rates and water production rates, respectively. Figure 3C shows the temperature distribution along the x-direction at y=77.5 m (see the dotted white line in Figure 3) on the 1800th production day. Figure 3D presents the average temperature of the global reservoir model. Figures 3E,F compare the pressure and temperature distribution between LGR and EDFM at 60-day production respectively. As one can see from Figure 3, the results of the proposed model agree well with the results of CMG STARS.

In addition to the single fracture case, the validation was conducted on a complex fracture network model (see Figure 4). As shown in Figure 4, there are 10 natural fractures as well as a hydraulic fracture in the hydrate reservoir model. The conductivity of the natural fractures used both in CMG and EDFM is 1,000 mD∙m.

Figure 5 compares the results obtained from EDFM against the results from CMG STARS, including gas production rate (Figure 5A), water production rate (Figure 5B), temperature distribution at y=77.5 m on the 1800th production day (Figure 5C), and the average temperature of the global reservoir model (Figure 5D). In addition, Figure 5e∼5F shows a comparison between LGR and EDFM on pressure and temperature distribution at 60-day production respectively. As shown in Figure 5, the results of the proposed model show slight differences from the results of CMG STARS.

Based on the results shown in Figure 3 and Figure 5, one can find that the proposed EDFM model is reliable in simulating the production of gas-hydrate reservoirs considering complex fracture networks and heat transfer. Thus, in this work, the authors will carry out a comprehensive study of the well production of fracture-filled hydrate reservoirs on the basis of this proposed EDFM model.

In order to explore the effect of fracture networks on the performance of gas-hydrate production well, three scenarios were investigated: Scenario #1 which considers 30 natural fractures and 1 hydraulic fracture; Scenario #2 which only considers 1 hydraulic fracture; and Scenario #3 which contains no fracture. Figure 6 shows the top view of three scenarios in this section. It is worth noting that, there are 2 natural fractures (green line in Figure 6A) connected with hydraulic fracture (red line in Figure 6A) for scenario#1. The hydraulic fracture length is 50 m, and the conductivity is 1,000 mD∙m for natural fractures and 2000 mD∙m for hydraulic fracture. The BHP is set to be a constant of 5 MPa.

Figure 7A compares the gas production rates of the three scenarios. In Figure 7A, one can find that the daily gas production rate (q _g ) and cumulative gas production (Q _g ) with natural and hydraulic fractures (scenario #1) are both higher than those only with one hydraulic fracture (scenario #2) and those without fractures (scenario #3). Comparing the cumulative production of scenario #2 to that of scenario #3, one can find that the hydraulic fracturing treatment can significantly improve gas well productivity. In addition, the existence of the natural fracture increases the cumulative production from 4.39 × 10⁶ m³ (scenario #2) to 6.39 × 10⁶ m³ (scenario #1). This indicates that both hydraulic fracture and natural fractures play important roles in influencing the performance of gas-hydrate wells.

Figure 7B shows simulation results of different scenarios at the end of the simulation. As shown in pressure map, a higher pressure drop can be observed along the hydraulic fracture and along the connected natural fractures (see scenario#1). This implies that natural fractures can noticeably increase the expanding speed of the drainage area of the production well due to its connection to hydraulic fractures, thus, leading to the highest well productivity among the three studied scenarios. Besides, a heterogeneous distribution can be observed in temperature field maps affected by natural and hydraulic fractures. The low-temperature zone gradually spreads over time along the fractures for the scenario with natural and hydraulic fractures (scenarios #1 and #2) while it gradually spreads along the radial direction for the scenario without fractures (scenario #3). In addition, the existence of natural and hydraulic fractures increases the hydrate dissociation rate. The maximum hydrate dissociation rate occurs near the hydraulic fracture, connected natural fractures, and the position of the wellbore.

In this section, different numbers of natural fractures (N _nf =2, 10, 30, and 50) and different numbers of connected natural fractures (N _{f_con} =0, 1, 2, 4) were examined to explore their influence on pressure and temperature distribution. As shown in Figure 8A, 0, 8, 28, and 48 natural fractures (blue lines in Figure 8A) are randomly generated in the hydrate reservoir, respectively. In addition to the randomly generated fractures, in all four cases, two natural fractures are deliberately generated to intersect with the hydraulic fracture. In addition, we further investigated the well productivity with different numbers of connected natural fractures (N _{f_con} ) to explore the effect of connected natural fractures on gas well productivity. The connected natural fractures indicate a natural fracture that is connected to the hydraulic fracture. As shown in Figure 8B, in the four studied reservoir models, the distribution and number of the randomly generated unconnected natural fracture (blue lines) remain unchanged through the four studied cases, whereas the number of the connected natural fracture (green lines) is increased from 0 to 4.

Table 2 summarizes the cumulative gas production of the different scenarios. As we can see in Table 2, the cumulative gas production rates are slightly increased as the N _nf is varied from 10 to 50. For example, the cumulative gas production is increased by 6.48% as N _nf is increased from 2 to 50. However, the cumulative gas production rate is significantly increased as the value of N _{f_con} is increased. In addition, the well productivity is more likely to be improved with less connected natural fractures. For example, the cumulative gas production is increased by 95.76 × 10⁴ m³ as N _{f_con} is increased from 0 to 1, whereas the increase is only 16.84 × 10⁴ m³ as N _{f_con} is increased from 2 to 4.

A comparison of pressure, temperature, and hydrate saturation after 1800 days of production with different values of N _nf and N _{f_con} is shown in Figure 9. It can be observed in Figure 9A that, although the natural fractures are increased from 2 to 50, the pressure, temperature, and hydrate saturation fields illustrate a negligible difference. This denotes that the natural fractures, which are not connected to the hydraulic fracture, exert a small influence on the productivity of the gas-hydrate well. In Figure 9B, one can find that the connected natural fractures can induce remarkable expansions of the drainage area. The expansion can be more significant as N _{f_con} is increased from 0 to 1 than that as N _{f_con} is increased from 2 to 4. The results in Table 2 and Figure 9 imply that the connected natural fractures improve the well productivity by increasing the drainage areas of the production well. The calculated results represent that the connected natural fractures can exert a more important influence on the well productivity than the unconnected natural fractures.

In this section, the authors discussed the impacts of dimensionless hydraulic fracture conductivity (C _fD ) on the productivity of a fractured vertical gas well. The dimensionless hydraulic fracture conductivity can be described as (Cinco et al., 1978): C f D = k f w f k m x f (32)

In Eq. 32, k _f and k _m is hydraulic fracture permeability and matrix permeability, respectively; x _f and w _f is hydraulic fracture half length and width, respectively.

Figure 10A displays the gas production rates of a fractured vertical well that is calculated with the proposed numerical model with different dimensionless hydraulic fracture conductivity. As we can see from Figure 10A, A larger dimensionless fracture conductivity leads to higher well productivity, which is attributed to a lower flow resistance in the fracture. One also can find that the q _g and Q _g express a slight difference if the dimensionless fracture conductivity exceeds 100. This is because the fracture can be regarded as infinite-conductive if the dimensionless fracture conductivity exceeds 100, and the flow resistance in the fracture can be ignored.

Figure 10B compares the field maps at the end of the 1800th day with different values of C _fD . As we can see from pressure maps in Figure 10B, the maximum pressure drop occurs at the center of the gas hydrate reservoir. As the fracture conductivity is increased, the expanding speed of the drainage area becomes higher (e.g., the yellow area of C _fD = 100 is larger than that of C _fD = 1 at the end of the 1800th day in pressure field maps). In addition, for a higher conductivity fracture, the temperature drops more rapidly near the fractures (see temperature field maps). This is attributed to the fact that the temperature is mainly influenced by thermal convection and solid-hydrate dissociation. Since the gas and water dissociated from solid hydrate is more readily to be produced through the fracture with a higher conductivity, more solid hydrate will be dissociated near the fracture. The solid-hydrate dissociation is endothermic which can result in a reduction of temperature. In hydrate saturation field maps, one can find that the hydrate saturation near the fracture is lower with higher fracture conductivity, which indicates that more solid-hydrate is dissociated near the fracture.

Moreover, the authors varied the conductivity of natural fractures together with the conductivity of the hydraulic fracture to compare their effects on well productivity. Table 3 summarizes the cumulative gas production for different hydraulic fracture conductivity and different natural fracture conductivity at the end of the simulation. In Table 3, both the natural fracture conductivity and the hydraulic fracture conductivity are varied from 25 mD∙m to 250 mD∙m. The calculated results show that the cumulative gas production is increased both as the hydraulic fracture conductivity and the natural fracture conductivity are increased. Besides, hydraulic fracture conductivity exerts a more significant influence on gas well productivity than natural fracture conductivity. For example, the maximum increase of the cumulative gas production that is induced by varying the natural fracture conductivity occurs with C _hf = 250 mD∙m, and the relative increase is 15.14% which is from 506.07 × 10⁴ m³ with C _nf =25 mD∙m to 582.69 × 10⁴ m³ with C _nf =250 mD∙m. However, with C _nf = 250 mD∙m, the relative increase of the cumulative production by varying the hydraulic fracture conductivity is 58.15% which is from 368.43 × 10⁴ m³ with C _hf =25 mD∙m to 582.69 × 10⁴ m³ with C _hf =250 mD∙m.

Hydrate dissociation rate also imposes a significant influence on the mass and heat transfer in developing gas hydrate. The kinetic model of Kim et al. showed that the hydrate decomposition rate is related to decomposition rate constant (K _D ), hydrate surface area (A _dec ), and pressure difference which can be expressed as, m h = K D A d e c P e − P g (33)

The decomposition rate constant K _D is defined as, K D = K d ⁡ exp − E R T (34) where K _d represents the intrinsic hydrate dissociation rate constant, mole/(day·KPa·m²). As we can see from Eqs 33, 34, the hydrate dissociation rate is a time-varying parameter in developing gas-hydrate. To quantitatively examine the effect of hydrate dissociation rate on the gas well productivity of fracture-filled hydrate reservoir, the different values of intrinsic hydrate dissociation rate constant (K _d ) were considered, including 0, 1.07 × 10⁴, 1.07 × 10⁹, and 1.07 × 10¹³ mol/(day·KPa·m²). It is noted that an intrinsic hydrate dissociation rate of 0 indicates that the hydrate dissociation is ignored, and the initial hydrate decomposition is more rapid as the increase of K _d with the same values of hydrate surface area and pressure difference. The gas production rates with different values of K _d are illustrated in Figure 11A. As we can see from Figure 11A, the plots show that the gas well productivity is increased as the hydrate dissociation rate is increased, indicating that hydrate dissociation is favorable for the production of gas-hydrate reservoirs.

Figure 11B shows the field maps with different values of K _d . An interesting observation is that the pressure drop is most noticeable with the smallest value of K _d = 0, which implies that the drainage area expands most rapidly without hydrate dissociation rate, whereas the temperature is reduced most significantly with the largest value of K _d = 1.07 × 10¹³ mol/(day·KPa·m²) (see temperature field maps in Figure 11B). It is worth noting that, in Sections 4.1–4.4, the highest well productivity comes together with the largest rate of drainage area expansion and the largest rate of temperature reduction. The calculated results in Figure 11B denote that the temperature change can be more highly related to the hydrate dissociation than the drainage area expansion. In addition, the hydrate dissociation rate can lead to a reduction of solid hydrate, thus, one can see that the hydrate saturation is lower near the fracture with a higher hydrate dissociation rate. The minimum matrix temperature along the natural fracture is around 6.3°C with the value of K _d = 1.07 × 10¹³ mol/(day·KPa·m²), whereas the temperature is 7.1°C with the value of K _d = 1.07 × 10⁴ mol/(day·KPa·m²). Furthermore, the thermotactic habit of hydrate nucleation may play a role in the secondary hydrate formation process which would be important in temperature distribution (Yang et al., 2023).