A new study of multi-phase mass and heat transfer in natural gas hydrate reservoir with an embedded discrete fracture model

Studies of the hydrate cores have shown that natural fractures can be frequently observed in hydrate reservoirs, resulting in a fracture-filled hydrate. Therefore, it is highly necessary for industries to predict the gas well productivity of fracture-filled hydrate reservoirs. In this work, an embedded discrete fracture model is applied to characterize the natural fractures of fracture-filled gas-hydrate reservoirs. The non-linear mass and energy conservation equations which are discretized with the finite-difference method are solved by the fully implicit approach, and the proposed model is justified by a commercial simulator. On the basis of the proposed model, we investigate the influences of natural fractures, fracture conductivity, and hydrate dissociation rate on the gas well productivity and the distributions of pressure, temperature, and hydrate saturation. The simulation results show that hydraulic and natural fractures exert significant impacts on the gas well productivity of the fracture-filled hydrate reservoirs, and the cumulative gas production is increased by 45.6% due to the existence of the connected natural fractures. The connected natural fractures can impose a more important influence on the gas well productivity than the unconnected natural fractures. The cumulative gas production is increased by 6.48% as N nf is increased from 2 to 50, whereas the increase is 43.38% as N f_con is increased from 0 to 4. In addition, A higher hydraulic fracture conductivity can be more favorable than a higher natural fracture conductivity for improving the gas well productivity, and a higher hydrate dissociation rate can lead to a lower temperature along fractures due to a more noticeable reduction of solid hydrate. This study provides a theoretical basis for developing fracture-filled hydrate reservoirs efficiently in the future.


Introduction
Natural gas hydrate (NGH) is formed from gas and water under particular pressure and temperature (Sloan et al., 2007), and 150-180 m 3 of gas can be approximately released from one volume of NGH at standard conditions (Carroll, 2020). At present, four different classes of gas hydrate reservoirs have been found in nature (Moridis, 2003;Moridis, 2004;Moridis et al., 2005). Class 1-2 hydrate reservoirs are composed of the overburden layer, gas hydrate layer, gaseous hydrocarbon or water layer, and under burden layer. Class 3 hydrate reservoir distinguishes from Classes 1 and 2 hydrate reservoirs due to the absence of the underlying mobile fluid layer. Class 4 hydrate reservoir is characterized by a low hydrate saturation (<10%) and disperse hydrate distribution (Moridis, 2003;Moridis, 2004;Moridis et al., 2005). Abundant natural gas hydrate resources have been observed in ocean margins and permafrost regions, and the estimated global technically recoverable gas hydrate is approximately 10 15 -10 18 m 3 (Carroll, 2020;Li et al., 2021b). At present, the depressurization method has been widely applied for developing gas hydrate (Liu et al., 2017). Numerous studies have been conducted to evaluate the potential of developing hydrate reservoirs ascribing to the fact that gas hydrates are known as one of the most promising clean energies (Boswell et al., 2012).
Developing gas hydrate is complicated which involves multiphase flow (gas-water) and hydrate phase transition. Moreover, thermal convection and hydrate decomposition occur during the development of gas-hydrate which can significantly affect temperature distribution and gas well productivity. Therefore, a comprehensive insight into the effect of hydrate phase transition on multi-phase mass and heat transfer of gas-hydrate reservoirs is one of the key issues for optimizing hydrate development. The kinetic model proposed by Kim et al. (1987) was widely used to calculate the hydrate decomposition rate caused by depressurization. In their model, the hydrate decomposition rate is related to hydrate surface area, gas fugacity, and temperature; Pooladi-Darvish and Hong. (2004)compared the effect of thermal conduction and thermal convection on the productivity of gas hydrate wells. Their simulation results showed that thermal conduction has a more significant influence in affecting heat transfer than thermal convection; Sun et al. (2005); Esmaeilzadeh et al. (2008) both established a numerical model considering the one-dimensional fluid flow to evaluate gas well productivity in developing gashydrate reservoirs by depressurization method, and the results of their studies showed fluid enthalpy and heat transfer coefficient are the most sensitive parameters in developing gas hydrate. Moridis, (2003), Moridis, (2004), Moridis et al. (2005). developed a new numerical submodule named Tough+Hydtate to study the mass and heat transfer of different types of hydrate reservoirs excluding Class 4. They claimed that Class 1 and Class 2 hydrate reservoirs hold significant promise as potential energy sources. Tonnet and Herri (2009) investigated the gas production and dissociation kinetics of hydrate using experimental and numerical methods. Their studies demonstrated that the mass transfers of gas and water impose a remarkable influence on hydrate decomposition for low permeability gas hydrate reservoirs. Zhao et al. (2014) developed a two-dimensional numerical model to analyze the effects of sensible heat transfer on gas well productivity of hydrate reservoirs. The calculated results of their works showed that the hydrate will dissociate rapidly at a high value of thermal conductivity in the early developing stage.
In all the aforementioned studies, the gas-hydrate reservoirs are assumed as a mono-porous media where the gas hydrate can only be observed in the pore space. However, the most recent studies show that natural fractures commonly exist in gas-hydrate reservoirs (Wang et al., 2011;Wang et al., 2019;Ning et al., 2020). Since natural fractures can significantly influence the gas well productivity of hydrate reservoirs, it is highly necessary for industries to develop a mathematical model to predict the gas well performance of the fracture-filled gas hydrate reservoirs considering the effect of natural fractures.
As one of the most applicable and economical method, numerical simulation has been widely used for evaluating the effect of the complex natural fracture network on gas well productivity. However, in terms of capturing the geometries of natural fractures, the traditional numerical simulator requires applying the local grid refinement (LGR) technique or unstructured grids, which can be inconvenient. In order to overcome the deficiencies of the traditional numerical methods, the embedded discrete fracture model (EDFM) method has been introduced to describe complex fractures (Lee et al., 2001). Compared to the methods of LGR and unstructured grids, the EDFM method characterizes the complex natural fractures with unrefined and structured grids. In recent studies, numerous EDFMs have been proposed to investigate the mechanism of mass and heat transfer in various naturally fractured reservoirs.
AlTwaijri et al., 2018 explicitly modeled complex fracture geometries by use of the EDFM method, and their simulation results confirmed that the complexity of the fracture network imposes a considerable influence in estimating gas and water recoveries. Li et al. (2021) combined an extended finite element method (XFEM) and EDFM to model the geothermal systems. The simulation results justified the application of EDFM for modeling fluid mass and heat transfer. Sangnimnuan et al., 2021 upgraded their numerical model by considering the geomechanics effects to model water-oil two-phase flow based on the EDFM. The calculated results indicated that the stress distribution and mass transfer are mainly affected by hydraulic fracture geometries ;Shi, (2021) combined the EDFM method with the boundary element method (BEM) to investigate the temperature change considering water-oil two-phase flow. Their calculated results showed that it is highly necessary to consider the temperature field effect on developing fractured geothermal system.
Although the EDFM was used widely to characterize the natural fracture in conventional oil/gas reservoirs, it has not been hitherto applied in gas-hydrate reservoirs to predict the gas well performance (Yan et al., 2016;Yan et al., 2018;Xu et al., 2017;Xu et al., 2018;Xu et al., 2021;Xu et al., 2022). For the scenarios of gas-hydrate reservoirs, multiphase flow, heat transfer, and hydrate dissociation all raise stringent barriers for one to build the EDFM model to evaluate the gas well productivity of the gashydrate reservoirs. Therefore, in this work, taking the effect of the natural fractures into consideration, we develop an EDFM model to characterize the gas well performance of the fracturefilled hydrate reservoirs. This proposed model is validated against a commercial simulator. Furthermore, with the aid of the proposed model, we analyze the influences of the complex fracture network, fracture conductivity, number of connected natural fractures, and hydrate dissociation rate on the gas well productivity and the distributions of pressure, temperature, and hydrate saturation.

Model assumptions
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; (2) the compressibility of rock and fluids are considered to characterize porosity change; (3) the reservoir is homogeneous and the flow of water and gas conforms to Darcy's law; (4) the capillary pressure is neglected in the model; (5) the effects of hydrate saturation changes, hydrate occurrence state, and pore structure on the relative permeability of gas and water are neglected; (6) only thermal conduction and thermal convection in the reservoir are considered, and the influence of thermal radiation is neglected.

Mass conservation equations
For the matrix system, the mass conservation equation of the gas and water can be expressed as (Hong and Pooladi-Darvish, 2003): 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,mf + ρ l q l,ff + m l z zt ρ l V v S l f , l g, w where ] l is volumetric flow rates, which can be expressed as: For hydrate, the mass conservation equations of the matrix and fracture systems are (Hong and Pooladi-Darvish, 2003): 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: 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: Likewise, the discrete form of fluid flow between fracture grids in Eq. 2 can be written as: The discrete form of mass transfer between matrix grid and fracture grid in Eqs 1, 2 can be written as: The discrete form of mass transfer between intersecting fracture grids in Eq. 2 can be written as: The discrete form of mass transfer between the wellbore and fracture grid in Eq. 2 can be written as: 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 Frontiers in Earth Science frontiersin.org #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: where T m k m ΔyΔz Δx 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): 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): 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): 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): For the hydrate, we can have:

Energy conservation equations
For the matrix system, the energy conservation equation can be expressed as (Hong and Pooladi-Darvish, 2003): In addition, the energy conservation equation for the fracture system is: 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

FIGURE 1
Schematic of the EDFM model for fractured hydrate reservoirs.
Frontiers in Earth Science frontiersin.org 04 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): 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), 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): 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): 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): Energy exchange U fw between well and the ith fracture can be expressed as: In general, the discrete form for the energy conservation equation of matrix and fracture grids can be expressed as (Song and Liang, 2009): 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).

Model validation
In order to examine the accuracy of the proposed model, we validate the proposed model against a commercial simulator CMG STARS. As the pressure is decreased in the reservoirs, the gas and water can be gradually released from the hydrate. To illustrate the fluid flow of different phases, we utilize the Brooks-Corey equation (Brooks, 1964) to characterize the relative permeability of water and gas which is expressed as: In this work, k rw,max =0.9, k rg,max =0.7, S wc =0.2, S gc =0.05, n w =n g =2. The other benchmark values of the parameters used in the simulation are listed in Table 1.

Hydrate reservoir with single fracture
In this case, we established a 155×155×30 m 3 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 Frontiers in Earth Science frontiersin.org 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.

Hydrate reservoir with complex fracture network
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 2
Top view of the grid system used in CMG STARS.
Frontiers in Earth Science frontiersin.org ( 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.

Results and discussion
With the aid of the proposed EDFM model, the authors studied the effects of natural fractures and heat transfer on the production of gas-hydrate reservoirs. The dimension of the reservoir is 310 × 310 ×   Table 1.

Effect of complex fracture network
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  Frontiers in Earth Science frontiersin.org

FIGURE 7
Simulation outputs of different scenarios: (A) Comparison of the production rates (B) pressure, temperature, and hydrate saturation field maps at the 1800th production day. Frontiers in Earth Science frontiersin.org 09 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 6 m 3 (scenario #2) to 6.39 × 10 6 m 3 (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.

Effect of number of natural fractures
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 4 m 3 as N f_con is increased from 0 to 1, whereas the increase is only 16.84 × 10 4 m 3 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.

Effect of fracture conductivity
In this section, the authors discussed the impacts of dimensionless hydraulic fracture conductivity (C fD ) on the productivity of a fractured

FIGURE 10
Simulation outputs of different scenarios with different values of hydraulic fracture conductivity: (A) Comparison of the production rates (B) pressure, temperature, and hydrate saturation field maps at the 1800th production day.
Frontiers in Earth Science frontiersin.org 11 vertical gas well. The dimensionless hydraulic fracture conductivity can be described as (Cinco et al., 1978): 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 infiniteconductive 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 4 m 3 with C nf =25 mD·m to 582.69 × 10 4 m 3 with C nf =250 mD·m.

FIGURE 11
Simulation outputs of different scenarios with different values of K d : (A) comparison of the production rates (B) pressure, temperature, and hydrate saturation field maps at the 1800th production day.
Frontiers in Earth Science frontiersin.org 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 4 m 3 with C hf =25 mD·m to 582.69 × 10 4 m 3 with C hf =250 mD·m.

Effect of hydrate dissociation rate
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, The decomposition rate constant K D is defined as, where K d represents the intrinsic hydrate dissociation rate constant, mole/(day·KPa·m 2 ). 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 4 , 1.07 × 10 9 , and 1.07 × 10 13 mol/(day·KPa·m 2 ). 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 13 mol/ (day·KPa·m 2 ) (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 13 mol/(day·KPa·m 2 ), whereas the temperature is 7.1°C with the value of K d = 1.07 × 10 4 mol/ (day·KPa·m 2 ). 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).

Conclusion
In this work, we developed an embedded discrete fracture model to study the gas well productivity in fracture-filled gas hydrate reservoirs. By using a finite-difference scheme and a fully implicit method to solve the non-linear equations, we validated the proposed model against a commercial simulator. In addition, we investigated the influences of natural fractures, fracture conductivity, and hydrate dissociation rate on mass and heat transfer of fracturefilled gas hydrate reservoir. The calculated results lead us to draw the following conclusions.
(1) This proposed numerical model is reliable for evaluating the gas well productivity of fracture-filled hydrate reservoirs considering the effect of natural fractures and heat transfer, and the research results are of great significance to guide the efficient development of fracture-filled hydrate reservoirs.
(2) Both hydraulic fractures and natural fractures can increase gas well productivity. Moreover, the connected natural fractures can exert a more important influence on the well productivity than the unconnected natural fractures. The cumulative gas production is increased by 6.48% as N nf is increased from 2 to 50, whereas the increase is 43.38% as N f_con is increased from 0 to 4. (3) The gas well production rate is increased by increasing fracture conductivity, and the conductivity of hydraulic fracture exerts a more significant influence on the gas well productivity than that of natural fracture. The relative cumulative gas production increment is 15.14% and 58.15% as C hf and C nf are increased from 25 mD to 250 mD, respectively. (4) The well productivity is increased as the hydrate dissociation rate is increased. Besides, a higher hydrate dissociation rate can lead to a lower temperature along fractures due to a larger reduction of solid hydrate. The minimum temperature along natural fractures is around 6.3°C with the value of K d = 1.07 × 10 13 mol/(day·KPa·m 2 ).

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Frontiers in Earth Science frontiersin.org