Multi-Lateral Well Productivity Evaluation Based on Three-Dimensional Heterogeneous Model in Nankai Trough, Japan

: Widely employed in hydrate exploitation, the single well method is utilized to broaden the scope of hydrate decomposition. Optimizing the well structure and production strategy is necessary to enhance gas recovery efﬁciency. Complex wells represented by the multilateral wells have great application potential in hydrate mining. This study focused on the impact of multilateral well production methods on productivity, taking the Nankai Trough in Japan as the study area. The spatial distribution of physical parameters such as porosity, permeability, and hydrate saturation in the Nankai Trough has signiﬁcant heterogeneity. For model accuracy, the Sklearn machine learning and Kriging interpolation methods were used to construct a three-dimensional heterogeneous geological model to describe the structure and physical property parameters in the study area of the hydrate reservoir. The numerical simulation model was solved using the TOUGH + Hydrate program and ﬁtted with the measured data of the trial production project to verify its reliability. Finally, we set the multilateral wells for hydrate high saturation area to predict the gas and water production of hydrate reservoir with different exploitation schemes. The main conclusions are as follows: 1 (cid:13) The Sklearn machine learning and Kriging interpolation methods can be used to construct a three-dimensional heterogeneous geological model for limited site data, and the ﬁtting effect of the heterogeneous numerical simulation model is better than that of the homogeneous numerical simulation model. 2 (cid:13) The multilateral well method can effectively increase the gas production rate from the hydrate reservoir compared with the traditional single well method by approximately 8000 m 3 /day on average (approximately 51.8%). 3 (cid:13) In the high saturation area, the number of branches of the multilateral well were set to 2, 3, and 4, and the gas production rate was increased by approximately 51.8%, 52.5%, and 53.5%. Considering economic consumption, the number of branching wells should be set at 2–3 in the same layer.


Introduction
Under natural conditions, methane hydrates are found in permafrost on continental shelves and in seafloor sediments [1,2]. There is about 160 volumes of methane for each volume of water at standard air pressure and temperature [3,4]. Therefore, the development and utilization of MH resources has become an attractive way for countries to deal with the energy crisis.
In sediments, solid hydrates are decomposed into liquid gas and water to produce hydrates. The decomposition of hydrates is achieved by raising the temperature above the equilibrium conditions, lowering the pressure below the equilibrium conditions, and changing the equilibrium conditions themselves. Makogon describes three common methods of hydrate separation: decompression, thermal injection, and suppression injection, which are based on traditional crude oil recovery techniques; decompression is the most common [5]. Theoretically, the in situ hydrate decomposition rate can be improved by increasing the bottom hole pressure drop. However, the strength of a hydrate reservoir is relatively low; blindly increasing the production pressure difference cannot raise the productivity effectively, but may cause overall collapse of the wellbore or reservoir failure resulting in sand burial of the borehole [6]. Therefore, enlarging the decomposition area is an effective alternative for increasing the gas production capability and alleviating the geotechnical problems caused by an excessive pressure drop.
Due to the large contact area between the multilateral well and the reservoir, the production can be effectively increased compared with the single well under the same pressure. Multilateral wells are produced from reservoir production [7,8]. Over the years, the number of multilateral completions has grown considerably due to the continued development of orientation and completion techniques [9]. So far, more than 8000 wells have been drilled worldwide, with remarkable economic benefits [10]. This technology has laid the foundation for the development and utilization of multi-well technology to develop gas hydrate.
Recently, the multilateral wells have been introduced into gas hydrate production wildly. Wilson et al. discovered in 2010 that, by adding multiple branches to a single vertical well, gas production in the North Slope of Alaska's gas hydrate reserves could be significantly improved [11], a finding that has been further supported by the recent introduction of multilateral wells into gas hydrate production. In the same year, Yamakawa et al. [12] and Wilson et al. [13] showed the benefits of a greater quantity of multilateral well branches in natural gas reservoir production from sandy reservoirs. Li et al. [14] proposed a technique with a large borehole and multi-lateral branches for the clayey silt hydrate reservoir in the South China Sea in 2019, with the primary goal of reducing sand output to maintain the reservoir's stability and significantly increase gas productivity. Lastly, Ning et al. [15] proposed multilateral well schemes with and without reservoir reformation technologies, such as hydraulic slotting and burden sealing, and the results demonstrate that the burden sealing is better for long-term gas recovery, while the multilateral well with hydraulic slotting is more suitable for short-term field trials. Li et al. proposed that complex structure wells, such as extended-reach horizontal and multilateral wells, should be taken into account when researching hydrate exploitation [16]. Moridis et al. [17] found that horizontal wells had a great advantage over vertical wells for gas production from Class 2 and Class 3 MH reservoirs. Furthermore, the hydrate in the mixed layer is always decomposed before that in the pure hydrate layer under any exploitation method, due to changes in hydrate solid concentration [18].
The above investigations conducted have been successful in advancing the utilization of multilateral wells in the production of natural gas hydrate. However, application of multilateral well in heterogeneous reservoirs were not discussed much. A more thorough exploration should be conducted into the most suitable amalgamation of geometric parameters, gas, and water production performance of a multilateral well in heterogeneous reservoir gas hydrate production.
The Japanese government has commenced a series of scientific investigations and probes into the creation of natural gas reservoirs. At Nankai Trough, the first offshore MH production test was conducted in March 2013 [19,20], but the gas production process only lasted six days and was brought to a close due to the intense sand production. Subsequently, the second offshore MH production test was conducted using depressurization in May to June of 2017 [21]. The resources of Nankai Trough require economical extraction and utilization, making the gas recovery enhancement an urgent issue. Despite the field tests' current gas production rates being far from the commercial exploration level of 5.0 × 10 6 m 3 /day [22], a multilateral well should be considered in this area to increase gas production in such complex geological conditions. At Site AT1 in the Nankai Trough, we conducted a study to explore gas production performance. Utilizing borehole geophysical logging, seismic interpretation, and core analyses, a geological model was created to compare a single vertical well under depressurization. Additionally, multilateral well models were constructed based on the single vertical well to explore the influence of geometric parameters and production interval manners on the production performance. It can be of great use for future gas production.

Characterization of Heterogeneous Reservoir Structure
The structure and physical parameters of the MH in the Nankai Trough exhibit a distinct spatial heterogeneity based on the data of previous projects. Identification of the hydrate reservoir is a crucial procedure for accurate simulation. However, traditional methods cannot easily characterize such an invisible underground space with complex structure. Research into the utilization of computer technology for numerical modeling is a relatively recent area of inquiry. In this study, machine learning is applied to hydrate reservoir identification. The multilayer perceptron (MLP) algorithm is a powerful algorithm in machine learning that can be used in classification problems. The Kriging interpolation method was used to accurately obtain the porosity, permeability, hydrate saturation, and physical parameters to construct a 3D heterogeneous geological model of the hydrate reservoir.

Geological Setting
At the Daini Atsumi Knoll off the coasts of the Atsumi and Shima peninsulas, Japan (Figure 1), the AT1 test site was located in the eastern Nankai Trough identified as one of the 16 methane hydrate concentrated zones (MHCZ). Figure 2 displays the depth profiles of lithofacies, hydrate saturation, porosity, and absolute permeability. Logging data and core sample analysis from the field test confirmed the presence of a 62 m thick MH reservoir at a depth of 1274-1336 m below sea level (mbsl), with an ocean depth of around 998 m [23]. The lithology of the MH reservoir was not homogeneous, but varied with depth, with three sub-zones of varying lithofacies being identified: upper sand/silt alternate layers (1274-1288 mbsl), middle silt-dominated layers (1288-1303 mbsl), and lower sand-dominated layers (1303-1336 mbsl) [24]. analyses, a geological model was created to compare a single vertical well under depressurization. Additionally, multilateral well models were constructed based on the single vertical well to explore the influence of geometric parameters and production interval manners on the production performance. It can be of great use for future gas production.

Characterization of Heterogeneous Reservoir Structure
The structure and physical parameters of the MH in the Nankai Trough exhibit a distinct spatial heterogeneity based on the data of previous projects. Identification of the hydrate reservoir is a crucial procedure for accurate simulation. However, traditional methods cannot easily characterize such an invisible underground space with complex structure. Research into the utilization of computer technology for numerical modeling is a relatively recent area of inquiry. In this study, machine learning is applied to hydrate reservoir identification. The multilayer perceptron (MLP) algorithm is a powerful algorithm in machine learning that can be used in classification problems. The Kriging interpolation method was used to accurately obtain the porosity, permeability, hydrate saturation, and physical parameters to construct a 3D heterogeneous geological model of the hydrate reservoir.

Geological Setting
At the Daini Atsumi Knoll off the coasts of the Atsumi and Shima peninsulas, Japan (Figure 1), the AT1 test site was located in the eastern Nankai Trough identified as one of the 16 methane hydrate concentrated zones (MHCZ). Figure 2 displays the depth profiles of lithofacies, hydrate saturation, porosity, and absolute permeability. Logging data and core sample analysis from the field test confirmed the presence of a 62 m thick MH reservoir at a depth of 1274-1336 m below sea level (mbsl), with an ocean depth of around 998 m [23]. The lithology of the MH reservoir was not homogeneous, but varied with depth, with three sub-zones of varying lithofacies being identified: upper sand/silt alternate layers (1274-1288 mbsl), middle silt-dominated layers (1288-1303 mbsl), and lower sanddominated layers (1303-1336 mbsl) [24].

Spatial Distribution of Hydrate Reservoir
Based on the distribution of several exploration and test production boreholes in the study area, the characterization range is the Daini Atsumi hill, where station AT1 is located.
The first step was to collect the drilling and seismic exploration data in the study area and obtain the coordinates, lithology, water depth, and other drilling data. The graphic processing modules VTK and Pyvista in Python were used to draw the distribution cloud map of the drilling and then draw the shape of the seafloor based on the distribution cloud map of the drilling hole ( Figure 3).  Based on the distribution of several exploration and test production boreholes in the study area, the characterization range is the Daini Atsumi hill, where station AT1 is located.
The first step was to collect the drilling and seismic exploration data in the study area and obtain the coordinates, lithology, water depth, and other drilling data. The graphic processing modules VTK and Pyvista in Python were used to draw the distribution cloud map of the drilling and then draw the shape of the seafloor based on the distribution cloud map of the drilling hole ( Figure 3).
In the second step, the Sklearn machine learning methods were used to characterize the spatial distribution of hydrate reservoirs. During characterization, the hydrate section information was obtained from the drilling hole. Log data were used to estimate the hydrate saturation out of the drilling hole area; the Archie equation used for estimation is (1).
R W is the resistivity of formation water; and S h is hydrate saturation; ϕ is density porosity, and R t is resistivity of actual formation. The main process is as follows: 1 Data preprocessing. The dataset was established from drilling, logging, and seismic exploration data, and the data were converted into matrix form by inputting the well location and lithology data. 2 Neural network. The MLP algorithm of Sklearn was used to train the model, and the hyperbolic tangent function (tanh) was selected as the activation function. The Determine the confusion matrix. We adopted the confusion_matrix function to determine the confusion matrix, which combines the predictions with actual conditions. 4 Set the size of the study area and conduct grid subdivision. The size of the study area was 2700 m × 2700 m × 150 m. A total of 40,000 grids were divided in the model, and the hydrate reservoir was identified. 5 Generate regular VTK grids and visualize the results.
The hydrate reservoir model established by the Sklearn machine learning method is shown in Figure 4a, and the red part represents the hydrate distribution. Figure 4a illustrates that the spatial distribution of the hydrate reservoirs is heterogeneous and discontinuous. Figure 4b is the seismic exploration data of the same section. In Figure 4b, the distribution shape of the actual MH is similar to the shape of seafloor. By comparing the two figures, the Sklearn machine learning method can be used to accurately simulate the shape of the MH, building a geological model that is similar to the actual structure of the MH.

Spatial Distribution of Hydrate Reservoir
Based on the distribution of several exploration and test pro study area, the characterization range is the Daini Atsumi hill, wh The first step was to collect the drilling and seismic explorat and obtain the coordinates, lithology, water depth, and other d processing modules VTK and Pyvista in Python were used to dra map of the drilling and then draw the shape of the seafloor based map of the drilling hole ( Figure 3).

Spatial Interpolation Methods
From the borehole and logging information, the physical property parameter distribution of multiple boreholes can be ascertained. Six surfaces of hydrate reservoirs below the seafloor were selected: 260, 278, 285, 291, 304, and 310 mbsf. Based on the spatial distribution obtained by machine learning, Kriging interpolation was used to perform porosity interpolation.
Kriging interpolation was employed to acquire the spatial distribution of hydrate saturation, permeability, and porosity in the vertical profile direction (2700 m × 1000 m) of the hydrate reservoir ( Figure 5). In Figure 5, there is a hydrate reservoir area with high saturation from 1300 to 2200 m in the X direction, and the hydrate reservoir thickness in this area is relatively large.

Spatial Interpolation Methods
From the borehole and logging information, the physical property parameter distribution of multiple boreholes can be ascertained. Six surfaces of hydrate reservoirs below the seafloor were selected: 260, 278, 285, 291, 304, and 310 mbsf. Based on the spatial distribution obtained by machine learning, Kriging interpolation was used to perform porosity interpolation.
Kriging interpolation was employed to acquire the spatial distribution of hydrate saturation, permeability, and porosity in the vertical profile direction (2700 m × 1000 m) of the hydrate reservoir ( Figure 5). In Figure 5, there is a hydrate reservoir area with high saturation from 1300 to 2200 m in the X direction, and the hydrate reservoir thickness in this area is relatively large.
The Sklearn machine learning and Kriging interpolation methods were used to obtain the spatial distribution of the hydrate reservoir structure, hydrate saturation, porosity, and permeability in the study area, which is more accurate than the homogeneous hydrate reservoir. The numerical simulation model can better reflect the production of gas and water behavior in the mining process, similar to the actual site. The model provides the basis for the prediction of depressurized mining productivity by multilateral wells.

Model Geometry and Spatial Discretization
Employing a cube of 900 m in length and breadth, with a thickness of 122 m for the model domain (Figure 3), 60 m of MHCZ (60 grids in Z direction) was used, as well as both overburden and underburden layers. A production well with a radius of 0.1 m was situated in the center of the model. Figure 6 displays a 38 m production interval from the top of MHCZ for the cross-section profile.
As for the mesh, based on previous numerical simulations, exploitation saw the principal processes occurring within a restricted area of the production [27]. For this reason, the grid density around the well should be increased appropriately. In the XY direction, the grid was dissected every 5 m within 100 m of the production well. The grid size was greater farther from the production well. Every 3 m in the MHCZ, it was dissected in the Z direction, and 2 m at the top of the MHCZ. Every 5 m, it was dissected in the overburden and underburden layers. The simulation domain was divided into 32 × 32 × 31 grids (31,744 cells) in the X, Y, and Z directions. The Sklearn machine learning and Kriging interpolation methods were used to obtain the spatial distribution of the hydrate reservoir structure, hydrate saturation, porosity, and permeability in the study area, which is more accurate than the homogeneous hydrate reservoir. The numerical simulation model can better reflect the production of gas and water behavior in the mining process, similar to the actual site. The model provides the basis for the prediction of depressurized mining productivity by multilateral wells. tain the spatial distribution of the hydrate reservoir structure, hydrate saturation, porosity, and permeability in the study area, which is more accurate than the homogeneous hydrate reservoir. The numerical simulation model can better reflect the production of gas and water behavior in the mining process, similar to the actual site. The model provides the basis for the prediction of depressurized mining productivity by multilateral wells.

Model Geometry and Spatial Discretization
Employing a cube of 900 m in length and breadth, with a thickness of 122 m for the model domain (Figure 3), 60 m of MHCZ (60 grids in Z direction) was used, as well as both overburden and underburden layers. A production well with a radius of 0.1 m was situated in the center of the model. Figure 6 displays a 38 m production interval from the top of MHCZ for the cross-section profile.
As for the mesh, based on previous numerical simulations, exploitation saw the principal processes occurring within a restricted area of the production [27]. For this reason, the grid density around the well should be increased appropriately. In the XY direction, the grid was dissected every 5 m within 100 m of the production well. The grid size was greater farther from the production well. Every 3 m in the MHCZ, it was dissected in the Z direction, and 2 m at the top of the MHCZ. Every 5 m, it was dissected in the overburden and underburden layers. The simulation domain was divided into 32 × 32 × 31 grids (31,744 cells) in the X, Y, and Z directions.  As for the mesh, based on previous numerical simulations, exploitation saw the principal processes occurring within a restricted area of the production [27]. For this reason, the grid density around the well should be increased appropriately. In the XY direction, the grid was dissected every 5 m within 100 m of the production well. The grid size was greater farther from the production well. Every 3 m in the MHCZ, it was dissected in the Z direction, and 2 m at the top of the MHCZ. Every 5 m, it was dissected in the overburden and underburden layers. The simulation domain was divided into 32 × 32 × 31 grids (31,744 cells) in the X, Y, and Z directions.
To compare with the heterogeneous model, a conceptual model in the homogeneous condition of the hydrate reservoir was established during the simulation. In the homogeneous condition, the conceptual model is an axisymmetric two-dimensional radial RZ2D model. The model's radius is 1000 m, and its Z-directional setting is identical to that of the heterogeneous model ( Figure 7). The production well is situated in the middle of the model, with its completion section atop the methane hydrate reservoir, measuring 38 m. In the radial direction, the production well is centered, and the farther from the production well, the larger the grid is. The minimum grid of 10 m around the production well was 0.5 m. Longitudinally, every 0.5 m, 244 grids were split, resulting in a total of 24,400 grids for the entire model. To compare with the heterogeneous model, a conceptual mod condition of the hydrate reservoir was established during the sim neous condition, the conceptual model is an axisymmetric two-di model. The model's radius is 1000 m, and its Z-directional setting i heterogeneous model (Figure 7). The production well is situate model, with its completion section atop the methane hydrate rese In the radial direction, the production well is centered, and the fa tion well, the larger the grid is. The minimum grid of 10 m arou was 0.5 m. Longitudinally, every 0.5 m, 244 grids were split, resu grids for the entire model.  Table 1 displays the main parameters of the heterogeneous termined by the spatial distribution of the reservoir physical pr dicted in the preceding section. This study employed the well-logg ple-analysis results from the field-trail set in Eastern Nankai Trou [28].
Under homogeneous conditions, the permeability of the form and hydrate saturation in the model have a homogeneous distribu   Table 1 displays the main parameters of the heterogeneous model, which were determined by the spatial distribution of the reservoir physical property parameters predicted in the preceding section. This study employed the well-logging data and core-sample-analysis results from the field-trail set in Eastern Nankai Trough of Japan extensively [28].

Reservoir Properties and Parameters
Under homogeneous conditions, the permeability of the formation, initial porosity, and hydrate saturation in the model have a homogeneous distribution (Table 2).

Initial and Boundary Conditions
The model varied in its pressure reduction, with a seafloor temperature of 3.75 • C. Before depressurization, the domain was in an equilibrium state. The lateral boundary 450 m from the production well was assumed to be a zero flow boundary, and both the overburden and underburden layers of the MHCZ hydrate reservoir had a boundary with constant pressure and temperature, which showed fluid and heat flow.
The Kriging interpolation results from the preceding section demonstrate that the porosity of the same layer remains relatively unchanged. Therefore, in the numerical model, the porosity distribution is considered as layered heterogeneity, and the pores of different layers are considered. Hydrate saturation refers to the hydrate reservoir distribution structure in the previous section and takes the average value of the corresponding thickness based on the hydrate saturation interpolation results. Figure 8 shows the gas hydrate saturation, initial porosity, and permeability. model, the porosity distribution is considered as layered heterogeneity, and the pores of different layers are considered. Hydrate saturation refers to the hydrate reservoir distribution structure in the previous section and takes the average value of the corresponding thickness based on the hydrate saturation interpolation results. Figure 8 shows the gas hydrate saturation, initial porosity, and permeability.

Comparison of Historical Fitting between Homogeneous and Heterogeneous Conditions
A comparison between the numerical simulation gas production rate under homogeneous and heterogeneous conditions and the measured data from the first engineering test production site of the Nankai Trough in Japan is depicted in Figure 9. Figure 9 shows the 6-day average gas production rates calculated by the heterogeneous and homogeneous models are 0.95 and 1.35 times that measured by the test project, respectively. During this period, there was rapid pressure reduction in the wellbore. After 3 days of mining, the disturbance effect on the hydrate reservoir around the wellbore during drilling decreased. The hydrate reservoir's decomposition during depressurization is the source of gas production at this juncture. Therefore, the numerical simulation results after 3 days of mining under heterogeneous conditions fit well with the actual site test results.

Comparison of Historical Fitting between Homogeneous and Heterogeneous Conditions
A comparison between the numerical simulation gas production rate under homogeneous and heterogeneous conditions and the measured data from the first engineering test production site of the Nankai Trough in Japan is depicted in Figure 9. Figure 9 shows the 6-day average gas production rates calculated by the heterogeneous and homogeneous models are 0.95 and 1.35 times that measured by the test project, respectively. During this period, there was rapid pressure reduction in the wellbore. After 3 days of mining, the disturbance effect on the hydrate reservoir around the wellbore during drilling decreased. The hydrate reservoir's decomposition during depressurization is the source of gas production at this juncture. Therefore, the numerical simulation results after 3 days of mining under heterogeneous conditions fit well with the actual site test results.

The Numerical Simulation Code
Using the TOUGH + HYDRATE simulator, a code for simulating hydrate bearing geologic systems, numerical simulations were conducted. This code solves coupled equations of mass and heat balance, allowing it to model the non-isothermal gas release, phase behavior, and flow of fluids and heat in complex geological media, from laboratory to reservoir, at any scale where Darcy's law is valid, in the typical natural methane hydrate deposits (e.g., permafrost and deep ocean sediments). bility (c).

Comparison of Historical Fitting between Homogeneous and Heterogeneous Conditions
A comparison between the numerical simulation gas production rate under homo geneous and heterogeneous conditions and the measured data from the first engineering test production site of the Nankai Trough in Japan is depicted in Figure 9. Figure 9 shows the 6-day average gas production rates calculated by the heterogeneous and homogene ous models are 0.95 and 1.35 times that measured by the test project, respectively. During this period, there was rapid pressure reduction in the wellbore. After 3 days of mining the disturbance effect on the hydrate reservoir around the wellbore during drilling decreased. The hydrate reservoir's decomposition during depressurization is the source o gas production at this juncture. Therefore, the numerical simulation results after 3 days of mining under heterogeneous conditions fit well with the actual site test results. Figure 9. Comparison of simulated gas production rate and measured data of test production engi neering under homogeneous and heterogeneous conditions.

The Numerical Simulation Code
Using the TOUGH + HYDRATE simulator, a code for simulating hydrate bearing geologic systems, numerical simulations were conducted. This code solves coupled equa tions of mass and heat balance, allowing it to model the non-isothermal gas release, phase

Design of a Single Vertical Well
In the center of the model, a 0.1 m radius vertical well was situated, with 38 m production interval from the MHCZ's summit. The wellbore pressure reduced rapidly from 10.3 to 4.4 MPa in the first 12 h. Then, the wellbore pressure slowly decreased before stabilizing at 4 MPa approximately 4.5 days after exploitation.

Design of Multilateral Wells
In this regard, we established several comparison cases. Based on the reports, when the branch angle was greater than 45 • , the productivity growth declined [31]. The maximum vertical depth of the lateral branch was limited when the dip angle of the lateral branch remained constant, taking into account the varying spatial distribution of hydrate saturation. The difference in gas production in MHCZ for limited thickness was minor when the difference in vertical depth of the lateral branch was a few meters; therefore, the effect of this parameter on the gas production effect of multilateral well was not considered in this study. At 45, we simulated the gas production performance of a multilateral well under depressurization.
In Section 2, Sklearn delineates the distribution of hydrate saturation in the formation. In this section, the influence of hydrate saturation on multilateral wells in heterogeneous layers is analyzed by changing the location of multilateral wells. The schemes are shown in Table 3. The single vertical well model mentioned previously was the basis for all parameters, including the initial and boundary conditions of multilateral vertical well models.  Figure 10 displays the development of the volumetric rate of gas and water production at a single vertical well, as well as the volumetric rates of water production over a year of production. Initially, the wellbore decompressed quickly and did not reach equilibrium (before t = 10 days). The gas production rate then surged and fluctuated greatly. From t = 10 days to t = 35 days, due to the widening of the decompression range, the gas hydrate sediment began to break down extensively, thus raising gas and water production. After 35 days of depressurization, the gas hydrate sediment decomposed as the temperature decreased, and the silt layer with a low permeability hindered the transmission of pressure and temperature in the vertical direction, leading to a decrease in gas production. This pressure difference caused water to flow into the production well, with an average rate of 2400 m 3 /day during one year of depressurization and 12,000 m 3 /day for gas production.

Evolution of Physical Properties
In Figure 11a, the physical properties in the hydrate reservoir, after a day of depressurization, 100 days, and 1 year, are displayed in a spatial distribution. The pressure decrease at the early stage was mainly concentrated around the production wells. The MHCZ2 layer was a hydrate reservoir dominated by silt with low permeability, and the low permeability layer hindered the rate of the pressure reduction in the Y and X directions. After one year of production, the range of depressurization (greater than 5 MPa) was approximately 30 m around the production well. The hydrate reservoir's permeability was greater, yet the transverse propagation range of pressure deduce was wider, and After 35 days of depressurization, the gas hydrate sediment decomposed as the temperature decreased, and the silt layer with a low permeability hindered the transmission of pressure and temperature in the vertical direction, leading to a decrease in gas production. This pressure difference caused water to flow into the production well, with an average rate of 2400 m 3 /day during one year of depressurization and 12,000 m 3 /day for gas production.

Evolution of Physical Properties
In Figure 11a, the physical properties in the hydrate reservoir, after a day of depressurization, 100 days, and 1 year, are displayed in a spatial distribution. The pressure decrease at the early stage was mainly concentrated around the production wells. The MHCZ2 layer was a hydrate reservoir dominated by silt with low permeability, and the low permeability layer hindered the rate of the pressure reduction in the Y and X directions. After one year of production, the range of depressurization (greater than 5 MPa) was approximately 30 m around the production well. The hydrate reservoir's permeability was greater, yet the transverse propagation range of pressure deduce was wider, and the effective pressure reduction (above 2 MPa) span was around 200 m. In contrast, in low permeability reservoirs, the transmission range of pressure reduction was small, the effective depressurization range was approximately 30 m, and the vertical transmission range of pressure reduction was blocked.
In Figure 11a, the physical properties in the hydrate reservoir, after a day of depres surization, 100 days, and 1 year, are displayed in a spatial distribution. The pressure de crease at the early stage was mainly concentrated around the production wells. Th MHCZ2 layer was a hydrate reservoir dominated by silt with low permeability, and th low permeability layer hindered the rate of the pressure reduction in the Y and X direc tions. After one year of production, the range of depressurization (greater than 5 MPa was approximately 30 m around the production well. The hydrate reservoir's permeabi ity was greater, yet the transverse propagation range of pressure deduce was wider, an the effective pressure reduction (above 2 MPa) span was around 200 m. In contrast, in low permeability reservoirs, the transmission range of pressure reduction was small, the e fective depressurization range was approximately 30 m, and the vertical transmissio range of pressure reduction was blocked. The thermal evolution of the hydrate reservoir, as depicted in Figure 11b, is cause by depressurization over the course of one day, a hundred days, and a year. Endotherm decomposition of the hydrate is the process that causes it to decay. As the hydrate d cayed, the temperature around the production well decreased in a similar fashion to th pressure, with a greater decrease in the high permeability area than in the low permeabi ity area. Figure 12 displays the one-year evolution of the volumetric rate (QG) of gas produce from various multilateral vertical wells compared to a single vertical well with depressu ization. It is evident that the multilateral well technique can significantly raise the ga 3 Figure 11. Evolution of physical properties in the hydrate reservoir. (a) Evolution of pressure in the hydrate reservoir due to depressurization for 1 day, 100 days, and 1 year. (b) Evolution of temperature in the hydrate reservoir due to depressurization for 1 day, 100 days, and 1 year.

Gas and Water Production Behaviors
The thermal evolution of the hydrate reservoir, as depicted in Figure 11b, is caused by depressurization over the course of one day, a hundred days, and a year. Endothermic decomposition of the hydrate is the process that causes it to decay. As the hydrate decayed, the temperature around the production well decreased in a similar fashion to the pressure, with a greater decrease in the high permeability area than in the low permeability area. surization. It is evident that the multilateral well technique can significantly raise the gas production rate from the hydrate reservoir by around 8000 m 3 /day on average, which is 51.8%, compared to the traditional single well method; however, no clear discrepancies were observed between the schemes of the two. The reasons for this probably as follows: the lower sand layer in the sand/silt alternate layers of the methane hydrate reservoirs may impede the downward transmission in pressure; thus, the vertical range of pressure and temperature in the multilateral well was hindered. Owing to the small difference in setting the lateral branch in vertical depth the under different schemes, the hydrate decomposition range of producing wells also varied less.

Gas and Water Production Behaviors
Energies 2023, 16, x FOR PEER REVIEW 14 of 1 production rate was approximately 800 m 3 /days higher for four branches on the same layer (Case6) than for two branches on the same floor (Case2). However, the average ga production rate was approximately 100 m 3 /day higher for four branches (Case6) than fo three branches (Case4) on the same layer; therefore, the higher number of branches did not improve the production. The graph in Figure 13 reveals a 1 year progression in the volumetric rate of wate production from various multilateral vertical wells compared to a single vertical well us ing depressurization. It is evident that the rate of water production during depressuriza tion is higher than that of a single vertical well, with an average of 1000 m 3 /day. This i due to the fact that the water produced was mainly created from the decomposition o hydrates in the hydrate reservoir during production; thus, the more hydrates decom posed, the more water was produced, and the water production rate was greater. The correlation between the water production rates is as follows: Case6 > Case4 > Case5 > Case2 > Case3 > Case1. The gas production rate of multilateral wells situated in the upper and lower layers of the hydrate reservoir (Case3 and Case5) was higher than those in the upper and middle layers of the hydrate reservoir with low permeability (Case4 and Case6). This was due to the high overall saturation and permeability of the upper layer (MHCZ1) and lower layer (MHCZ3), which enabled the transmission of temperature and pressure at the early stage and the decomposition of hydrate. However, the longer the distance between the two, the smaller the decrease in temperature and pressure, which was not conducive to hydrate exploitation in the later stage. From 10 days to 300 days, the same number of branches and similar vertical depth of the lateral branch were used to measure this. At the latter stage of extraction, the hydrate reservoir around the production well had finished decomposing, and the temperature and pressure decreased along X and Y directions due to the obstruction of the underlying sand layer. In the middle of the MHCZ2, the overall hydrate saturation and permeability were low. Therefore, this arrangement was beneficial to hydrate exploitation in the later stage. At the same location, with a similar vertical depth of the lateral branch, when there were more branches in the multilateral well, the depressurization superposition effect was improved, increasing the gas production rate. The gas production rate was approximately 800 m 3 /days higher for four branches on the same layer (Case6) than for two branches on the same floor (Case2). However, the average gas production rate was approximately 100 m 3 /day higher for four branches (Case6) than for three branches (Case4) on the same layer; therefore, the higher number of branches did not improve the production.
The graph in Figure 13 reveals a 1 year progression in the volumetric rate of water production from various multilateral vertical wells compared to a single vertical well using depressurization. It is evident that the rate of water production during depressurization is higher than that of a single vertical well, with an average of 1000 m 3 /day. This is due to the fact that the water produced was mainly created from the decomposition of hydrates in the hydrate reservoir during production; thus, the more hydrates decomposed, the more water was produced, and the water production rate was greater. The correlation between the water production rates is as follows: Case6 > Case4 > Case5 > Case2 > Case3 > Case1. Figure 12. Volution of volumetric gas production rates for different cases of multilateral wells during the 1 year production period.
The graph in Figure 13 reveals a 1 year progression in the volumetric rate of water production from various multilateral vertical wells compared to a single vertical well using depressurization. It is evident that the rate of water production during depressurization is higher than that of a single vertical well, with an average of 1000 m 3 /day. This is due to the fact that the water produced was mainly created from the decomposition of hydrates in the hydrate reservoir during production; thus, the more hydrates decomposed, the more water was produced, and the water production rate was greater. The correlation between the water production rates is as follows: Case6 > Case4 > Case5 > Case2 > Case3 > Case1. Figure 13. Evolution of volumetric water production rates for different cases of multilateral wells during the 1 year production period. Figure 13. Evolution of volumetric water production rates for different cases of multilateral wells during the 1 year production period.

Evolution of Physical Properties
The spatial distribution of hydrate saturation in the hydrate reservoir over the course of a year of depressurization production is depicted in Figure 14, with various multilateral well scenarios. Due to the low permeability layer located in the model, which hindered the vertical depressurization rate, the hydrate in MHCZ1 decomposed more than MHCZ3. Thus, the MHCZ3 decomposition of Case3 and Case5 is sufficient, which may be because multiple branching wells are distributed with high hydrate saturation. Through comparative analysis of Case1, Case3, and Case5, we found that when multiple branching wells are distributed in the same position, the higher the number of branching wells is, the more complete the decomposition will be; however, when the number of branching wells exceeds 4, the increase is not significant.

Comparative Analysis
By utilizing the multilateral well method, the gas production volume and rate from the hydrate reservoir can be significantly increased in comparison to the traditional single well approach. Furthermore, the multilateral branches based on the main wellbore can be augmented, thereby increasing the hydrate decomposition contact area, and thus raising the pressure and temperature transmission efficiency.
Compared with the base case, the gas production volume can be increased by an average of 41% when multilateral wells are located in a hydrate high-saturation area. In the high-saturation area, the number of branches of the multilateral well were set to 2, 3, and 4, and the gas production rate was increased by approximately 51.8%, 52.5%, and 53.5%. Thus, there is an increasing trend, but it is not obvious.
In Figure 15, the gas/water ratio (R GW ) of a single vertical well over a one-year production period is depicted, as well as the various cases of multilateral wells. R GW is a relative criterion for evaluating the hydrate production efficiency. The R GW after 1 year of production in a multilateral well was higher than that in a single vertical well. However, there was little difference in R GW between the different multilateral well cases. The inset shows a comparison of R GW for different multilateral well cases after 10 to 100 days of depressurization production. Case3 (three branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) and Case5 (four branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) had relatively large R GW during this period.
The spatial distribution of hydrate saturation in the hydrate reservoir over the cou of a year of depressurization production is depicted in Figure 14, with various multilate well scenarios. Due to the low permeability layer located in the model, which hinder the vertical depressurization rate, the hydrate in MHCZ1 decomposed more th MHCZ3. Thus, the MHCZ3 decomposition of Case3 and Case5 is sufficient, which m be because multiple branching wells are distributed with high hydrate saturatio Through comparative analysis of Case1, Case3, and Case5, we found that when multi branching wells are distributed in the same position, the higher the number of branchi wells is, the more complete the decomposition will be; however, when the number branching wells exceeds 4, the increase is not significant.

Comparative Analysis
By utilizing the multilateral well method, the gas production volume and rate fro the hydrate reservoir can be significantly increased in comparison to the traditional sin well approach. Furthermore, the multilateral branches based on the main wellbore can augmented, thereby increasing the hydrate decomposition contact area, and thus raisi the pressure and temperature transmission efficiency.
Compared with the base case, the gas production volume can be increased by an a erage of 41% when multilateral wells are located in a hydrate high-saturation area. In t high-saturation area, the number of branches of the multilateral well were set to 2, 3, a In the two cases with the highest gas-to-water ratio, a comparative analysis was conducted based on their gas production as shown in Figure 16. Case3 increased production by 40.9%, and Case5 increased production by 41.7% (Figure 16). tive criterion for evaluating the hydrate production efficiency. The RGW after 1 year of production in a multilateral well was higher than that in a single vertical well. However, there was little difference in RGW between the different multilateral well cases. The inset shows a comparison of RGW for different multilateral well cases after 10 to 100 days of depressurization production. Case3 (three branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) and Case5 (four branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) had relatively large RGW during this period. Figure 15. Evolution of gas-to-water ratio (RGW) of different cases of multilateral wells during the 1 year production period. The inset shows the evolution of RGW from 0 to 100 days.
In the two cases with the highest gas-to-water ratio, a comparative analysis was conducted based on their gas production as shown in Figure 16. Case3 increased production by 40.9%, and Case5 increased production by 41.7% (Figure 16).  ization production. Case3 (three branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) and Case5 (four branches per layer, branching in high-permeability reservoirs MHCZ1 and MHCZ3) had relatively large RGW during this period. Figure 15. Evolution of gas-to-water ratio (RGW) of different cases of multilateral wells during the 1 year production period. The inset shows the evolution of RGW from 0 to 100 days.
In the two cases with the highest gas-to-water ratio, a comparative analysis was conducted based on their gas production as shown in Figure 16. Case3 increased production by 40.9%, and Case5 increased production by 41.7% ( Figure 16).

Conclusions and Suggestions
We employed the Nankai Trough in Japan as our research site. Utilizing the geological information from Site AT1, a combination of well-logging and core-sample data was used to create a numerical simulation model for the extraction of hydrate reservoirs, which was then solved with the TOUGH + Hydrate program; the gas extraction from hydrate reservoirs under different extraction conditions was studied, and we predicted the gasproduction potential from the sedimentary-complex-hydrate reservoir. We can draw the following main conclusions.
(1) The Sklearn machine learning and Kriging interpolation methods can be used to construct a three-dimensional heterogeneous geological model for limited site data, and the fitting effect of the heterogeneous numerical simulation model is better than that of the homogeneous numerical simulation model. (2) In the same layer, the gas production rate increased by approximately 51.8%, 52.5%, and 53.5% when the number of branches in the multilateral well were 2, 3, and 4, respectively. Although there was an increase, it was not significant. Considering the economic cost, the number of branches of the multilateral well should be set at 2-3.
(3) The layout of multilateral wells is affected by reservoir heterogeneity. When the multilateral wells are located in the formation with high hydrate saturation, the hydrate decomposition is improved. (4) The optimization of the lateral geometric parameters proved intricate in the field test.
Additionally, it should take into account indicators such as drilling cost and estimated production duration. To further enhance gas production efficiency, we can investigate the effect of other parameters on the production of multilateral wells.