The Fetkovich analytic aquifer model is based on the principle of material balance; that is, the volume of water at the initial time under the current pressure should be equal to the sum of the pore volume under the current pressure and the volume of water invading the adjacent reservoir (Fetkovich, 1971). When the source and sink are not considered in the aquifer, they can be expressed as W B w − [ W B w i − c f ( p i − p ¯ ) W B w i ] = W e (1) Equation 1 can be transformed to: p ¯ = − ( p i W e i ) W e + p i (2) The rate of invasion water can be expressed as follows: q w = J w ( p ¯ − p w f ) (3) By combining Eq.2 and Eq.3, the following equation is obtained: W e = W e i p i ( p i − p w f ) { 1 − e − [ ( q w i ) max / W e i ] t } (4) According to the above equations, the total water invasion rate in the time interval Δ t n can be obtained as follows: Δ W e = W e i p i ( p ¯ a , n − 1 − p ¯ n ) { 1 − e − [ ( q w i ) max / W e i ] Δ t n } (5)

At present, the commonly used numerical simulation tools for HBL development include the STOMP-HYD simulator of the Pacific Northwest National Laboratory in the United States, the MH21-HYDRES simulator of Japan, and the TOUGH + HYDRATE and HydrateResSim simulators of the Berkeley National Laboratory in the United States. (Li et al., 2016; Luo et al., 2020). Among these simulators, the phase equilibrium and kinetic models are the most widely used models (Moridis, 2003; Wan et al., 2020). The equilibrium model considers the hydrate formation and dissociation to occur at chemical equilibrium, and the system is always assumed to be equilibrium. For laboratory scale models, the timestep is quite small (seconds or minutes) and the assumption of equilibrium is hard to be satisfied. Therefore, the kinetic model should be selected under such circumstances. However, the timestep of field scale model is days or even months. The assumption of equilibrium is easily satisfied, and therefore, the equilibrium and kinetic model obtain quite similar results (Kowalsky and Moridis, 2007). The calculation process of the kinetic model is more complex, so the convergence is relatively poor, and the calculation time is longer. Therefore, a phase equilibrium model was used in this study considering that the field scale model is used. The components used only include methane and water, and the mass conservation equation for each component can be expressed as follows: ∂ ∂ t ( ∑ β = w , g ϕ S β ρ β X β κ ) + ∇ ⋅ ( ∑ β = w , g X β κ F → β ) = q κ (6) The energy conservation equation is as follows: ∇ ⋅ { − [ ( 1 − ϕ ) K r + ∑ β = w , g ϕ S β K β ] ∇ T + ∑ β = w , g H β F → β + ε σ ∇ T 4 } + ∂ ∂ t [ ( 1 − ϕ ) ρ r C r T + ∑ β = w , g ϕ S β ρ β H β + H d ] = q h (7) Hydrate formation and dissociation is based on the phase equilibrium condition. The phase equilibrium curve used in this study is the regression form proposed by Moridis et al. (2011b) as follows: T ≥ 273.2 K         ln  P e = { - 1 . 94138504464560 × 10 5 + 3.31018213397926 × 10 3 T - 2 . 25540264493806 × 10 1 T 2 + 7.67559117787059 × 10 − 2 T 3 - 1 . 30465829788791 × 10 − 4 T 4 + 8.86065316687571 × 10 − 8 T 5 T < 273.2 K         ln  P e = { - 4 . 38921173434628 × 10 1 + 7.76302133739303 × 10 − 1 T - 7 . 27291427030502 × 10 − 3 T 2 + 3.85413985900724 × 10 − 5 T 3 - 1 . 03669656828834 × 10 − 7 T 4 + 1.09882180475307 × 10 − 10 T 5 (8)

Figure 1 shows a diagram of the coupling process of the Fetkovich aquifer model and the HBL model. As shown in this figure, the widely used two-way coupling method was used to couple the two models. First, the HBL and aquifer models were initialized, and the temperature of the aquifer was set to be the same as the initial temperature of the adjacent HBL grid. When the simulation began, the pressure, temperature, and hydrate saturation of the HBL were obtained according to the numerical simulation model of the HBL. When the pressure of the grids adjacent to the aquifer is lower than the initial pressure of the aquifer, the total water invasion rate from the aquifer to the HBL in the next time step is calculated using Eq. 5, and the aquifer pressure is updated according to Eq. 2. According to the principles of the numerical simulation, the mass flow rate of water from the aquifer to the grids adjacent to the aquifer can be calculated according to the following equation: q w j = ρ w Δ W e V B Δ t n W I D j ( λ w B w + λ g B g ) j ∑ j = 1 N [ W I D ( λ w B w + λ g B g ) ] j (9) W I D is given by (Moridis, 2003): W I D = 2 π k h ln r e r w + S − 3 4 (10) The heat of the invading water was obtained according to the water temperature and the calculated invasion rate. Then, the mass flow rate and heat of the invasion water are considered the source and sink terms of mass and energy, respectively, in Eq. 6 and Eq. 7 to calculate the pressure, temperature, and fluid saturation of the HBL in the next step. The above steps are repeated until the predetermined simulation time T _S is reached.

No similar coupling model has been reported in the literature; therefore, no simulation results can be directly used to validate the coupled model proposed in this article. However, in the two-way coupling process, only related data are transferred, and the procedure used to compute outputs from the two models is not changed (Tran et al., 2005). Therefore, if the water invasion rate is computed correctly, the results can be considered reliable. Many commercial software packages of the petroleum industry have incorporated the Fetkovich model. In this study, the commercial software CMG-STARS module was used to validate the calculation results of the coupled model. To make the simulation results comparable, a single water phase that can be simulated by both the STARS module and the coupled model is used. The coupled model built is shown in Figure 2A and its basic parameters are listed in Table 1. As shown in Figure 2A, the model uses radial coordinates and is divided into 28 vertical grids, of which the top three grids represent the overburden with a total thickness of 10 m, and the remaining 25 grids represent a water layer with a total thickness of 25 m. According to the basic parameters in Table 1, the pore volume of the water layer is 9.42 × 105 m³. The analytical aquifer is connected to the bottom of the water layer, and the outer boundaries are all no-flux boundaries except the bottom (water in the aquifer will flow from the aquifer to the water layer through the bottom). The initial pressures of the water layer and the analytical aquifer were both 13 MPa. A vertical well is located in the center and produces water at a fixed pressure of 4 MPa, which means that the inner boundary (wellbore boundary) is the Dirichlet boundary.

Three cases were compared in which the pore volumes of the aquifer were 4, 20, and 100 times that of the bottom water layer. Figure 2B shows a comparison of the water production rates obtained by the STARS module and the coupled model. It can be seen that the initial water production rate is almost the same for the three cases, which indicates that the elastic energy of the water layer is sufficient at the beginning. However, with an increase in time, the difference in water production rates among the different cases gradually increases. An aquifer with a larger pore volume can provide more water flow from the aquifer to the adjacent layer, and therefore, the larger the pore volume of the aquifer, the slower the deceleration of the water production rate. From the comparison, it can be seen that the water production rates of the coupled model established in this study are consistent with those calculated by the STARS module, and the coupled model is reliable.

The advantage of the coupled model is that it uses an analytical aquifer model to characterize the bottom water layer, while the single HBL model can only increase the pore volume of the bottom water layer by increasing the number of grids. Therefore, the calculation time of the single model increased significantly when the aquifer was large. Figure 3 shows a comparison of the calculation results and calculation time of the coupling model and the single model. Because the case in which the pore volume of the aquifer is 100 times that of the bottom water layer requires too many grids for the single model and the simulation time is too long, only the other two cases are compared. As shown in Figure 3A, the water production rates simulated by the coupled and single models are almost the same, which indicates that the analytical aquifer model can accurately simulate the water invasion, and the calculation results of the coupled model are reliable. As shown in Figure 3B, when the pore volume of the aquifer is four times that of the bottom water layer, the calculation time of the coupled model is approximately 2 min (CPU: Intel Core i7 7,700; Memory: 8 GB), while the calculation time of the single model is approximately 7 min. When the pore volume of the aquifer increases to 20 times that of the bottom water layer, the calculation time of the coupled model is approximately 3 min, which is not a significant increase. However, the calculation time of the single model was increased to 56 min owing to the increase in the grid number. The computation time of the single model is 18.7 times that of the coupled model. The above results show that the coupled model can ensure calculation accuracy and significantly reduce the calculation time.

The numerical simulation model was demonstrated using the basic parameters of the HBL at station SH7 in the Shenhu area of the South China Sea, and the initial temperature and pressure of the model were 14.15°C and 13 MPa, respectively (Liu et al., 2012; Li et al., 2011). According to Eq. 8, the equilibrium pressure corresponding to 14.15°C is 11.53 MPa, and therefore, the hydrate is stable under the initial conditions. The model can be divided into overburden, hydrate, and bottom water layers as shown in Figure 4. The thickness of the overburden layer is 10 m and that of the hydrate layer is 20 m. It is recommended that aquifers be modeled at least partially using water-filled grid blocks because if flow reversal occurs, the model with no water-filled grid blocks will encounter convergence problems [CMG (Computer Modelling Group), 2020]. The aquifer used in this study is composed of two parts. One part is the bottom water layer with 10 m thickness and the other part is an analytical aquifer represented by the Fetkovich model. A vertical well was located at the center of the formation. Similar to the model built in Validation and Computational Efficiency Analysis of the Coupled Model Section, the outer boundaries are all no-flux boundaries except the bottom, and the inner boundary is the Dirichlet boundary. The vertical well maintained a constant bottomhole pressure of 4 MPa during the entire depressurization process, and the perforated section of the vertical well ran through the entire hydrate layer. In the model, the hydrate phase cannot flow, and its relative permeability is 0. The modified STONE model is used for the gas-water relative permeability, which is expressed as follows (Moridis, 2003; Liu et al., 2020): k r W = min { [ S W − S i r W 1 − S i r W ] n W , 1 } k r G = min { [ S G − S i r G 1 − S i r A ] n G , 1 } (11) The widely used van Genuchten function is used for the capillary pressure, and it can be expressed as follows (Hou et al., 2016): P c a p = − P c o [ ( S W − S i r W 1 − S i r W ) − 1 λ − 1 ] 1 − λ (12) The specific parameters of the model are listed in Table 2. Three cases were simulated with the analytical aquifer model, and the total pore volumes of the aquifer (the sum of the pore volumes of the water layer and the aquifer model) were 4 (case 2), 20 (case 3), and 100 (case 4) times that of the hydrate layer. A class III HBL model (case 1) which only has overburden, underburden, and hydrate layers was also simulated for comparison. This case does not have a water layer or an analytic aquifer model; therefore, the pore volumes of the water layer and aquifer model are 0. The simulation time for all the four cases was 10 years.

The gas production rate and cumulative gas production for each case are shown in Figure 5. It can be seen that the pore volume of the aquifer has a significant impact on the depressurization performance. The peak daily gas production rates for cases 1 to 4 are 4.24 × 10⁴, 3.33 × 10⁴, 2.61 × 10⁴, and 1.20 × 104 m³/d, respectively, and the peak gas production times are 250, 310, 550, and 1,270 days, respectively, indicating that the higher the pore volume of the aquifer, the lower the gas production peak, and the later the gas production peak appears. This is mainly because when the aquifer is large, the amount of water that can flow into the HBL from the aquifer is also large, and the depressurization rate of the HBL is accordingly low. The hydrate dissociation rate is positively correlated with the amplitude of depressurization; therefore, the larger the aquifer, the lower the dissociation and gas production rates. In the later stage of depressurization development, the gas production rate shows the opposite trend; that is, the larger the aquifer, the higher the gas production rate. This is mainly because hydrate dissociation is an endothermic reaction. When there is no aquifer, hydrate dissociation in the early stage of depressurization causes the reservoir temperature to rapidly decrease and move closer to the phase equilibrium temperature of the hydrate; thus, the dissociation rate decreases significantly, and the gas production decreases accordingly. When the bottom water layer exists, the heat energy contained in the invasion water can promote hydrate dissociation, and therefore, the gas production for case 4 with the largest aquifer is also the highest in the later stages of depressurization. However, it can be seen from Figure 5B that case 1 without an aquifer can achieve the highest cumulative gas production in 10 years of depressurization development, followed by cases 2, 3, and 4 in order, where the largest aquifer achieves the lowest cumulative gas production.

A comparison of the water production rate and cumulative water production for each case is shown in Figure 6. It can be seen that in the early stage of depressurization, the water production rates for all the four cases were high. However, when there is no bottom water (case 1), the water production rate decreases rapidly with a decrease in the dissociation rate of the hydrate, while in the cases with bottom water, the decrease in the water production rate is slower. In the late stage of depressurization, the water production rates for cases 2 and 3 are close to those for case 1 owing to the limited pore volume of the aquifer. However, because of the large pore volume of the aquifer and sufficient water supply, the water production rate for case 4 is still significantly higher than those for the other cases after 10 years of depressurization development.

Figure 7 shows changes in the average pressure, and Figure 8 shows the pressure distribution of each case after 10 years of depressurization. It can be observed that the average pressure in all cases decreases with time. However, the larger the aquifer volume, the slower the pressure drop. After years of depressurization, the pressure of the entire HBL dropped to approximately 4MPa when there was no bottom water, as shown in Figure 8A. The pressure distributions for case 2 (Figure 8B) and case 3 (Figure 8C) are similar to those for case 1 owing to the smaller pore volume of the aquifer. However, when the pore volume of the aquifer was large (Figure 8D), only the reservoir pressure near the well was close to the bottom hole pressure, while the reservoir pressure far away from the well was still high, which indicates that the HBL can achieve effective depressurization when the aquifer is small, but it is difficult to achieve rapid depressurization when the bottom water is sufficient.

Figure 9 shows the changes in average temperature, and Figure 10 shows the temperature distribution for each case after 10 years of depressurization. It can be seen that the temperatures for all the cases decreased owing to the dissociation of the hydrate. The temperature of the whole HBL dropped to approximately 4.2°C, which corresponds to the phase equilibrium temperature of the bottom hole pressure of 4 MPa at the end of depressurization development because there is no external energy supplement for case 1 (Figure 10A). For cases 2 and 3 (Figures 10B,C), the temperature of the area far away from the vertical well also dropped to near the phase equilibrium temperature, but there was an evident cone-shaped high-temperature area at the bottom of the well. This is mainly because bottom water coning on the vertical well can supplement the heat energy consumed by hydrate dissociation near the bottom of the well. For case 4, the temperature of the HBL was significantly higher than the phase equilibrium temperature. This is mainly because it is difficult for the HBL to achieve rapid depressurization when the aquifer is large (Figure 10D); thus, the hydrate dissociation rate is lowered, and the thermal energy of the HBL is not fully utilized. Similar to cases 2 and 3, the temperature of the bottom water coning area in Figure 10D is also significantly higher than that of the non-coning area because of the energy supplement of the bottom water.

Figure 11 shows the changes in average hydrate saturation of the hydrate layer and Figure 12 shows the hydrate saturation distribution for each case after 10 years of depressurization. When depressurization begins, the hydrate near the well quickly dissociates because the near-well zone has the lowest pressure and the average hydrate saturation decreases. At the end of the depressurization development, the hydrate in the near-well zone was completely dissociated in all cases. However, the shapes of the hydrate dissociation region of the four cases were different. For case 1, the distribution of the hydrate saturation is approximately symmetrical in the vertical direction, while for other cases with bottom water, there is a cone-shaped dissociation zone near the bottom water; in addition, the larger the pore volume of the aquifer, the larger the area of the cone-shaped dissociation zone. This is mainly because the heat carried by the bottom water can effectively promote hydrate dissociation in the water invasion area. At the end of depressurization, there is still a large quantity of un-dissociated hydrate in each case, and the larger the aquifer, the greater the saturation of the remaining hydrate.

From the above analysis, it can be seen that bottom water can have a significant impact on the performance of the HBL. The larger the aquifer, the slower the depressurization rate, and the lower the corresponding gas production. Therefore, the invasion of bottom water is not conducive to the depressurization development of Class II HBLs.