Representative Dynamic Accumulation of Hydrate-Bearing Sediments in Gas Chimney System since 30 Kyr BP in the QiongDongNan Area, Northern South China Sea

: A stratigraphic complex composed of mass transport deposits (MTDs), where the gas occurrence allows for the formation of a gas chimney and pipe structure, is identified based on seismic interpretation in the QiongDongNan area of the northern South China Sea. During the Fifth Gas Hydrate Drilling Expedition of the Guangzhou Marine Geological Survey, this type of complex morphology that has close interaction with local gas hydrate (GH) distribution was eventually confirmed. A flow-reaction model is built to explore the spatial–temporal matching evolution process of massive GH reservoirs since 30 kyr before the present (BP). Five time snapshots, including 30, 20, 10, and 5 kyr BP, as well as the present, have been selected to exhibit key strata-evolving information. The results of in situ tensile estimation imply fracturing emergence occurs mostly at 5 kyr BP. Six other environmental scenarios and three cases of paleo-hydrate existence have been compared. The results almost coincide with field GH distribution below the bottom MTD from drilling reports, and state layer fracturing behaviors always feed and probably propagate in shallow sediments. It can be concluded that this complex system with 10% pre-existing hydrates results in the exact distribution and occurrence in local fine-grained silty clay layers adjacent to upper MTDs.


Introduction
Gas hydrate (GH) is an ice-like nonstoichiometric crystalline solid formed of water and methane.Natural GH can be found on the seabed, in ocean sediments, in deep lake sediments, as well as in the permafrost regions.The amount of carbon potentially trapped in natural hydrate-bearing deposits may be between 3 × 10 15 to 1 × 10 16 m 3 , which makes them of major interest as a potential energy resource [1,2].The presence and occurrence of GH in oceanic continental slope areas have attracted extensive attention and particular interest for the sake of energy potential, marine geohazards, and climate change issues [3,4].
A seismically imaged gas chimney is a type of stratigraphic sedimentary and hydraulic area where vertical seismic reflection profiles generally indicate that seafloor geological structures in this layer possess excellent transport capacity for upward gas-bearing fluids, and also lots of methane gas is trapped in conduit pores [5][6][7].It ordinarily extends several strata sequences, provides an effective aisle for functionally connecting deep hydrocarbon sources with the shallow bottom boundary of the local gas hydrate stability zone (HSZ), and conveys these multiphase fluids into the seafloor surface through affiliated pipe structures on many occasions [8,9].Hovland and Sommerville (1985) reported the identification of gas-charged sediments and the characteristics of two natural gas seepages in the North Sea, which suggested a gas-liquid mixture fluid seepage caused the pockmark development [10].After that, the impacts of seabed pockmarks associated with gas chimneys on offshore construction and gas hydrate occurrence have been extensively investigated by Hovland and collaborators [11,12].In passive continental marginal slope regions, diverse deposition activities along the down-slope direction promote shear stress to exceed the shear strength and subsequently cause sedimentary failure, thus generating another laterally extended complicated structure of mass transport deposits (MTDs) [13,14].
Plentiful biogenic methane produced from shallow organic matter and thermogenic methane produced by deep source rocks cross together through a main gas chimney area and enter into shallow seabed deposits, contributing to the formation of complex structural seepage GH reservoirs in a suitable thermodynamic environment [15][16][17].If the fluent pipe structure is disturbed or even stabbed across by MTDs, the upward flow of methanebearing fluids is blocked by a sudden decrease in layer permeability and is thus inclined to stop and remain in local layers to form abundant GHs with various occurrence states.Subsequently, pore pressure below gradually increases, and hydraulic failure is likely to occur.Therefore, gas chimneys and local MTD systems play a key role in the dynamic migration, occurrence, and accumulation of GH reservoirs [18][19][20].These activities are presumed to have occurred in many gas chimney-hydrate geological systems from time to time with geologic and environmental impacts [21].
Combining the visible distribution states of GH from the drilling core samples, we discovered that the accumulation type of GH reservoirs in the QiongDongNan Basin (QDNB) area of the northwestern South China Sea (SCS) (Figure 1a), to some extent, is similar to other typical marine reservoirs, such as the Ulleung basin in South Korea and the Krishna-Godavari (KG) basin in India [22,23].The GH system in the QiongDongNan (QDN) area is deemed to have the same distribution and accumulation traits described above.Previous gas and oil resource exploration surveys have identified large-scale gas chimney morphotypes in the target sediments by analyzing high-resolution seismic sections.The large range of low-frequency chaotic seismic reflection area, which is obviously different from the surrounding stratigraphy, indicates the presence of abundant gas.The boundary of this chimney scope was confined by the termination of the blanking reflection.The instantaneous frequency and some of the internal events can clearly reveal the fluid pathway of gas chimney and pipe structure.Bottom simulating reflector (BSR) characteristics usually appear at the top of these shapes, and fluid leakage channels show a pull-up trace of seismic indications reaching the seafloor [24][25][26].During the Fifth Gas Hydrate Drilling Expedition of the Guangzhou Marine Geological Survey (GMGS5) expedition, visible GH crystals recovered from drilling cores mainly show various features such as massive, laminates, nodules and veins massive, layers, nodules, and veins (sites W07, 08, and 09).GHs cement into fractured clay is the dominant phase in fine-grain silty sediments and have a concentrated distribution [27,28].In recent years, remarkable progress has been made in understanding how local seismic chimneys and MTDs control fluid migration and determine the growth patterns of GH reservoirs from geologic, geochemical, and geophysical studies according to field investigation and the drilling core analysis results of GMGS5 [29,30].In this study, we analyze GH reaction kinetics into the frame of a spatio-temporal dynamic evolution scale in In recent years, remarkable progress has been made in understanding how local seismic chimneys and MTDs control fluid migration and determine the growth patterns of GH reservoirs from geologic, geochemical, and geophysical studies according to field investigation and the drilling core analysis results of GMGS5 [29,30].In this study, we analyze GH reaction kinetics into the frame of a spatio-temporal dynamic evolution scale in a well constrained geological background.We considered that the evolution process of GH reservoirs is also controlled by environmental factors such as seafloor deposition, eustatic sea-level change, and organic gas supply.We use a flow-reaction model to explore how the interactions of sedimentary properties, fluid qualities, and reaction kinetics can lead to the current occurrence and distribution of GHs.We also focus on the growth and characteristics of fracture pathways in sediments within HSZ and how it propagates over GH content.

Geological Setting
The GMGS5 investigation area is located in the central channel of the upper slope of the QDN basin, adjacent to the Xisha Trough to the northeast (Figure 1b,c).Since the Neogene, the QDN basin has been in a post-rift thermal subsidence tectonic evolution stage.Wide-spreading Miocene (Huangliu formation), Pliocene (Yinggehai formation), and Quaternary (Ledong formation) marine sedimentary sequences formed in this period (Figure 1d).Since then, the local high deposition rate has led to the maximum sedimentary accumulation thickness along slope surface stretches about several kilometers [31][32][33][34].The detrital sediments of the Holocene and Pleistocene age have a sand content of more than 50%.Elastoplastic unconsolidated argillaceous particles fill into skeleton matrixes.They serve as effective barriers to stop fluid from further smoothing upward seepage and provide sufficient space for GH storage.Previous exploration revealed that deep hydrocarbon source rocks developed below the Huangliu formation and have entered the mature hydrocarbon threshold.The thermal evolution degree of organic matter is in the range of a mature and high mature light oil and gas window, which has the geological conditions of forming thermogenic gas.Moreover, the strata above the Huangliu formation are within the formation condition of biogenic gas because most of the source rocks developed in the low mature-immature stage [35][36][37][38].Intensive biogenic gas manifestation occur in petroleum drilling boreholes, especially those shallower than 2300 m below seafloor (mbsf).Both coalmeasure mudstone (which generates thermogenic methane) and marine mud shale (which produces biogenic methane extensively) are potential source rocks here.The measured carbon isotope of methane is −43.65-39.92‰ the mixture origin of gas source [30].
A large number of gas chimneys are found in the study area, with a prevalent northeast direction.The acoustic blanking area is mushroom-shaped with the appearance status of a small lower part and a large upper part [24,39,40].Most of these structures are firmly associated with the structures of faults and fractures [41].Almost all gas chimneys here vertically extended from a T3 to T1 sequence, and their tops can be up to 200-150 mbsf gauging from two-way travel time [36].Through characterizing the weak amplitude and chaotic seismic sections, three sets of MTD since the Late Quaternary have been confirmed above T1 interface, and laterally stretch across gas chimneys.With a high amplitude and reversed polarity, BSR here tends to sit above gas chimneys and just below the third MTD.Its wave impedance difference between the upper and lower layers is large and also shows a blank zone, indicating that probably abundant free gas is available within this study seabed [29,42].

GHs in Site W08
The 2D seismic profile and instantaneous frequency profile collected in the GMGS5 expedition have been carefully inspected.From the characteristics of the BSRs, logging while drillings (LWDs), and lithologic coring results, the site W08 can be deemed as a representative well of this complex system because it is located exactly above the gas chimney section (Figure 2a,b).The deep seismic chimney is estimated to originate from 3 Ma BP.It is elliptic plane-shaped and has a diameter of several kilometers.The upward hydrocarbon fluids are reckoned to pass through this section and reach the T1 interface.The top of this seismic chimney shows the BSR occurrence and nearly coincides with the bottom of MTD3.A shallow fluid pipe shape in the seismic profile is observed, which has a diameter of 180 m.It is columnar and calculated to be merely 0.5 Ma years old.This pipe is located just above the gas chimney area.Fluids will likely escape to the seafloor through it because the dome and pockmarks on the local seafloor that surface from remotely operated vehicle (ROV) observation indicate signs of fluid activity.Local high sedimentation rates since the early Pliocene easily trigger disequilibrium pore compression and bring about the overpressure outcome.In many situations, the fluid flow is mainly driven by the pore overpressure [24,40].The lithology of the MTD3 here has a 20-50% clay content which is inferred to be older than 0.5 Ma BP.GHs in drilling core images mainly appear messy and discontinuous, and they grow in the fractures within MTD3.With strata stress and thermodynamic condition change, the performance of gas-driven hydraulic fracturing is presumed to occur extensively in this GH system.The top of this seismic chimney shows the BSR occurrence and nearly coincides with the bottom of MTD3.A shallow fluid pipe shape in the seismic profile is observed, which has a diameter of 180 m.It is columnar and calculated to be merely 0.5 Ma years old.This pipe is located just above the gas chimney area.Fluids will likely escape to the seafloor through it because the dome and pockmarks on the local seafloor that surface from remotely operated vehicle (ROV) observation indicate signs of fluid activity.Local high sedimentation rates since the early Pliocene easily trigger disequilibrium pore compression and bring about the overpressure outcome.In many situations, the fluid flow is mainly driven by the pore overpressure [24,40].The lithology of the MTD3 here has a 20-50% clay content which is inferred to be older than 0.5 Ma BP.GHs in drilling core images mainly appear messy and discontinuous, and they grow in the fractures within MTD3.With strata stress and thermodynamic condition change, the performance of gas-driven hydraulic fracturing is presumed to occur extensively in this GH system.

Sea-Level Change and Deposition
The depth of seawater bottom in the GMGS5 investigation zone ranges from 1600 to 1800 m below sea level (mbsl).The deepest bottom boundary where GH core samples are drilled out is 237 mbsf.Nearly all GH layers occur above T1 interface (Figure 2b,c).The basin subsides since the Quaternary are roughly steady, and the sedimentary evolution is mainly controlled by the eustatic sea level change which is closely related to the Last Glacial Cycles that ended about 15,000 years ago.The detrital input is brought into the marginal sea by paleo-rivers [43].Since seabed position changes can be calculated according to relative sea level and deposition rate vsed (cm/kyr) (Figure 3a), the temperature and

Sea-Level Change and Deposition
The depth of seawater bottom in the GMGS5 investigation zone ranges from 1600 to 1800 m below sea level (mbsl).The deepest bottom boundary where GH core samples are drilled out is 237 mbsf.Nearly all GH layers occur above T1 interface (Figure 2b,c).The basin subsides since the Quaternary are roughly steady, and the sedimentary evolution is mainly controlled by the eustatic sea level change which is closely related to the Last Glacial Cycles that ended about 15,000 years ago.The detrital input is brought into the marginal sea by paleo-rivers [43].Since seabed position changes can be calculated according to relative sea level and deposition rate v sed (cm/kyr) (Figure 3a), the temperature and hydrostatic pressure at the bottom of the seawater are determined accordingly (Figure 3b).The recent rise in sea level in the area begins at about 30 kyr BP from historical inversion.hydrostatic pressure at the bottom of the seawater are determined accordingly (Figure 3b).The recent rise in sea level in the area begins at about 30 kyr BP from historical inversion.The distribution of deposit particle sizes and the skeleton in the seabed are often chaotic, and in particular, they suffer from the MTD influence.The measured local porosity distribution is plotted in Figure 3c.In order to simplify the treatment on the anomalistic sedimentary layers, initial porosity ϕ0 in deposits without GH existence is gauged by field data and its expression with depth z (mbsf): 0 ( ) 0.29 ( 0.29) exp( 0.00505 ) where porosity at the sedimentary surface ϕtop is 0.62.After GH forms in porous soil, we assume that the alteration of new porosity ϕh is suitable: where Sh is GH saturation.The distribution of deposit particle sizes and the skeleton in the seabed are often chaotic, and in particular, they suffer from the MTD influence.The measured local porosity distribution is plotted in Figure 3c.In order to simplify the treatment on the anomalistic sedimentary layers, initial porosity ϕ 0 in deposits without GH existence is gauged by field data and its expression with depth z (mbsf): where porosity at the sedimentary surface ϕ top is 0.62.After GH forms in porous soil, we assume that the alteration of new porosity ϕ h is suitable: where S h is GH saturation.

Methane in Pore Water
Through calculating the phase equilibrium curves of methane in liquid and gas states under local thermodynamic conditions [45][46][47], we compared how in situ GH is vertically distributed at GMGS5 site W08 (Figure 3d).It can be seen that methane supplies below HSZ bottom and even in the low part of HSZ are pretty abundant, so the current free gas is likely to be plentiful in these zones.We also inferred that maybe the GH reaction is still ongoing, while the methane gas in the middle and upper sections is completely converted into hydrates.The consequence is that merely GH and liquid water are preserved in local pores.According to the data from degassing cores [29], the initial methane concentration n diss (mmol/kg) in pores is quantified as: n diss = exp(4.17+ 0.007 × z − 5.18 × 10 −6 × z 2 ). (3)

Sedimentary Soil-Water Properties
We continue to adopt traditional normalized relating content to gauge sedimentary soil-water characteristics curves (SWCC), which is [48,49] as follows: where S e is effective fluid saturation; S l and S g are liquid and gas saturation, respectively; S li and S gi are initial irreducible liquid and gas saturation, and here S gi is always 0. Irreducible liquid water S li is considered to be controlled by clay content x c and ϕ h as [50]: The saturation-dependent relative permeability for gas k rg and liquid k rl phases can be defined by Corey's model: and where the pore size profile exponent π is 7.3 in this simulation [51].As the pressure difference between pore gas P g (MPa) and liquid phases P l (MPa), capillary pressure P c (kPa) is defined by pore entry pressure P th (kPa): where ω is another Brooks-Corey parameter.The demarcation on capillary drainage and the imbibition curves of shallow seabed sediments has been reasonably introduced by Daigle et al. (2020) [52] when tracing how gas-driven tensile fractures in marine hydratebearing layers.We use the same measurement to show the influence of clay fractions x c , which are: and The clay fraction scope in its seabed soil is 25-60%.When setting the layer porosity to be minimum 0.3 and maximum 0.7, the scope of the corresponding S li ranges from 0.2781 to 0.3788.

Hydraulic Conductivity in Gas Chimney
Compared with the sedimentary transportation quality of fluids in other parts of northern SCS [53], we can use an exponential expression to depict the strata intrinsic permeability k h (mD) after GH existence in the QDN area as: where k 0 (mD) is layer instinct permeability without GH existence.Based on the above description of the chimney-hydrate system, we set two different permeability values to reflect MTD influence.Within the 0-180 mbsf section, the permeability k 0 is reckoned as 2 mD.Below 180 mbsf, it is just within the gas chimney area and its value is 200 mD [30,36].

Tensile Strength Estimation
We also abide by the assumption that tensile failure can only occur when pore pressure steps across the resultant of minimum principle stress and the tensile strength of seabed soil together [52,54].The final force difference criterion term S v (MPa) in this work is: The vertical stress σ v (MPa) and horizontal stress σ h (MPa) can be linked according to Poisson's ratio in the QDN area, and the expression is: The sub-seafloor σ v is simplified as σ v = ρ b gz, and here the bulk sedimentary density ρ b is 1910 kg/m 3 .The tensile strength T sp (MPa) is determined as: where the Hoek-Brown constant m i is 4, and compressional wave velocity V p (m/s) at site W08 can be retrieved from the in situ logging curve as:

GH layer Formation Dynamics
From the comprehensive perspectives of gas source, kinetic reaction, fluid migration pathway, HSZ alteration, MTDs, and sedimentary soil properties, we established a onedimensional dynamic evolution model, which is in line with the previous geological background, to demonstrate the accumulation process of hydrate-bearing sediments and put specific parsing steps and corresponding expressions in the Appendix A of this paper.
When the seawater depth is larger than 500 mbsl, the relationship of temperature at the seafloor surface T sf ( • C) with the depth L w (mbsl) is [55]: The local temperature T se ( • C) in sediments is determined by geothermal gradient T d ( • C/m) and T sf ( • C), which is: All the parameter values in this basic simulation scene are listed in the following Table 1.
Table 1.Parameters used for deducing the basic scenario here.

Evolution Process
The old seafloor depth at 30 kyr BP is calculated at1790 mbsl in the light of the current sea level position.The earliest HSZ bottom boundary is 147 mbsf at that time.After the sedimentary process piles up, porewater pressure in accumulated sediments increases slightly to 0.0101 MPa/m, related to hydrostatic pressure 0.01 MP/m [56][57][58].To better exhibit how this complex system consisting of a gas chimney, GH, pipe, and MTDs evolves, the following six key controlling factors, including strata pressure, temperature, capillary pressure, phase saturations, permeability, and salinity were chosen to represent system characteristics, as they can provide enough insightful strata information.Table 1 tells the basic thermodynamic and hydraulic environment about this complex GH system.Both the initial gas and hydrate content are 0 in this basic scene.
Five moments including four divided times (30, 20, 10, and 5 kyr BP), and the present, are selected to show the system evolution at that time.In Figure 4a, the positive pressure change is set as the difference between sedimentary pore pressure and seawater hydrostatic pressure if there is no deposition activity.Similarly, in Figure 4b, positive temperature change is defined as the difference between sedimentary temperature and seawater temperature if there are no sediments.Both of them verify the previous judgment that the depth of the local seafloor position increases from 30 to 5 kyr BP and then decreases a bit until now.GH is calculated to accumulate considerably within the lower part of HSZ and the maximum saturation of GH from model simulation can reach 47% (Figure 4c).It can be seen that both the GH distribution close to the HSZ bottom area and its peak of saturation coincide with the drilling data very well.Simulation results reveal that free gas will gradually gather at the top zone of the gas chimney, which is also the bottom zone of HSZ and MTD3.The current gas saturation is 20.3%, to possibly provoke BSR signs.Our results also demonstrate that GH and free gas can hardly coexist in pores, according to the strata setting of the gas chimney quality.The alternation of salinity and capillary pressure has a strong positive correlation with GH saturation, while permeability displays a wholly negative correlation (Figure 4d-f).These phenomena are consistent with the general cognition of marine GH sediments.Overall, the flow-reaction process, strata features, and the combination of related key parameter values perfectly explain how gas chimneys and MTD determine the dynamic evolution and final occurrence of GH systems.

Force Disequilibrium
The distribution of sedimentary stress and strength and corresponding force balance in this geological system attracted widespread attention in the academic community.Whether strata hydraulic fracturing occurs in clay-dominated fine-grained layers and how cracks are created by pore overpressure propagating in coarse-grained turbidities are essential to understanding GH existence and occurrence.Based on the stress and strength

Force Disequilibrium
The distribution of sedimentary stress and strength and corresponding force balance in this geological system attracted widespread attention in the academic community.Whether strata hydraulic fracturing occurs in clay-dominated fine-grained layers and how cracks are created by pore overpressure propagating in coarse-grained turbidities are essential to understanding GH existence and occurrence.Based on the stress and strength analysis in Equations ( 12)-( 15), we calculated their values at these five moments in Figure 5a.The results show that the local tensile strength gradually decreases from the largest value at 30 kyr BP to the minimum at 5 kyr BP.After that, it increases rapidly.The current tensile strength is the same as that at 10 kyr BP.Meanwhile, the results of the criterion term, S v , indicate that except for the bottom area of the studied strata at the beginning time of 30 kyr BP, the phenomenon of tensile fracture will always occur (Figure 5b).The strata are most likely to break at 5 kyr BP, and the existence of GH can greatly improve the possibility of fracturing behavior in the strata.GH saturation has an observably positive change with the term S v .

Different Sedimentary Conditions
Local gas supply, kinetic reaction, and permeability are three decisive variables in our model to simulate this GH system.We designed six other geological scenarios that can represent possible combinations of strata and fluid environments to cultivate diverse occurrences of hydrate-bearing sediments (Table 2).We presented the differences between the results of free gas, GH, and fracturing criterion term in these six scenarios and the previous basic condition (Figure 6a,b).Scenario 1 represents extremely high permeability in the gas chimney and the upper strata.The results show that the GH distribution range is nearly the same, but the saturation peak can only reach 15%.The gas peak can be up to 85%, and its distribution range is wider.GH and free gas can coexist in the HSZ.Scenario 2 shows less gas supply and ki-

Different Sedimentary Conditions
Local gas supply, kinetic reaction, and permeability are three decisive variables in our model to simulate this GH system.We designed six other geological scenarios that can represent possible combinations of strata and fluid environments to cultivate diverse occurrences of hydrate-bearing sediments (Table 2).We presented the differences between the results of free gas, GH, and fracturing criterion term in these six scenarios and the previous basic condition (Figure 6a,b).
Table 2. Six possible geological scenarios for GH generation that may occur in the QDN area. in the HSZ with a peak value of 23%.Scenario 6 represents the largest gas supply, a smaller reaction coefficient, and a larger layer permeability in the gas chimney.GHs can fill the entire HSZ where sporadic free gas exists in the pores.By calculating these six fracturing criterion terms Sv, all the values are positive, indicating that fracturing will form or occur (Figure 6c).Meanwhile, the greater the hydrate content, the more likely that hydraulic fracturing will happen.

Palaeo-GH Existence
There is also a possibility that a certain amount of GH has already existed in the original HSZ before 30 kyr BP.With the rather big shift of the HSZ bottom boundary, the GH in the lower part decomposes until it vanishes, while the upper part remains there.We set the three initial hydrate contents as 5% (case1), 10% (case2), and 15% (case3) to investigate the evolution process since 30 kyr BP.Similarly, we employ the difference of GH, free gas and fracturing criterion term Sv between these three cases and the corresponding value in the basic case.From the results, it can be seen that free gas contents also increase accordingly (Figure 7a).GH peaks at the HSZ bottom in the three cases can reach 53%, 57% and 59%, respectively.The possibility of a second peak is 9%, 19%, and 31% at about 105 mbsf where they nearly coincide with MTD2 bottom in Figure 7b.The criterion term curves are still positive, indicating that strata fracture will develop in these settings (Figure 7c).Furthermore, the higher the hydrate content, the more intensive the fracture behavior.On the whole, it is believed that the initial GH saturation of 10% is nearly perfect and more in line with the current observation about the GH existence close to MTD2 and MTD3 from drilling logging.Scenario 1 represents extremely high permeability in the gas chimney and the upper strata.The results show that the GH distribution range is nearly the same, but the saturation peak can only reach 15%.The gas peak can be up to 85%, and its distribution range is wider.GH and free gas can coexist in the HSZ.Scenario 2 shows less gas supply and kinetic reaction coefficient.The content of hydrates and free gas both decline dramatically.Scenario 3 shows a larger gas supply and a smaller reaction coefficient.The results indicate that the distribution of GH is wider.The highest position can reach about 95 mbsf, but the GH peak which shifts upward can only reach 23%.The free gas content is the same, but the peak position is also shifted higher.Scenario 4 represents a higher gas supply and the slowest reaction coefficient, while the upper sedimentary permeability increases.The results show that although the hydrate and free gas are more widely distributed, the contents of them are very low.Scenario 5 represents a higher gas supply and a slower kinetic reaction coefficient.The results show that although the hydrate distribution range is wide and the peak is high enough, there is also a wide range of free gas content in the HSZ with a peak value of 23%.Scenario 6 represents the largest gas supply, a smaller reaction coefficient, and a larger layer permeability in the gas chimney.GHs can fill the entire HSZ where sporadic free gas exists in the pores.By calculating these six fracturing criterion terms S v , all the values are positive, indicating that fracturing will form or occur (Figure 6c).Meanwhile, the greater the hydrate content, the more likely that hydraulic fracturing will happen.

Palaeo-GH Existence
There is also a possibility that a certain amount of GH has already existed in the original HSZ before 30 kyr BP.With the rather big shift of the HSZ bottom boundary, the GH in the lower part decomposes until it vanishes, while the upper part remains there.
We set the three initial hydrate contents as 5% (case1), 10% (case2), and 15% (case3) to investigate the evolution process since 30 kyr BP.Similarly, we employ the difference of GH, free gas and fracturing criterion term S v between these three cases and the corresponding value in the basic case.From the results, it can be seen that free gas contents also increase accordingly (Figure 7a).GH peaks at the HSZ bottom in the three cases can reach 53%, 57% and 59%, respectively.The possibility of a second peak is 9%, 19%, and 31% at about 105 mbsf where they nearly coincide with MTD2 bottom in Figure 7b.The criterion term curves are still positive, indicating that strata fracture will develop in these settings (Figure 7c).Furthermore, the higher the hydrate content, the more intensive the fracture behavior.On the whole, it is believed that the initial GH saturation of 10% is nearly perfect and more in line with the current observation about the GH existence close to MTD2 and MTD3 from drilling logging.

Petrophysical Variations
Except for the heterogeneity of sedimentary layers, other lithological variations of mixtures in the marine subsurface, including pore shape and connection, particle size distribution, and clay contents, exert strong impacts on fluid transport ability and the stiffness of skeleton endmembers.In our model, although the uniform lithology assumption is used to gauge this deposition-flow-reaction process, different values of the intrinsic permeability in separate gas chimneys and pipe zones have been set to roughly reflect the real lithological situation.One of the two uncertainties lies in the usage of the empirical parameters of SWCC properties in Equations ( 5)-(11).In our simulation, the general values of sedimentary factors are adopted, which define the drainage-imbibition performance of seabed soils.Another is the failure criterion to judge distinct Mohr-Coulomb hydraulic fracturing.We do not elaborate on the scopes of Poisson's ratio, clay content, and the Hoek-Brown constant here, because they are the mean values and can explain most ordinary cases in our research deposits.

Petrophysical Variations
Except for the heterogeneity of sedimentary layers, other lithological variations of mixtures in the marine subsurface, including pore shape and connection, particle size distribution, and clay contents, exert strong impacts on fluid transport ability and the stiffness of skeleton endmembers.In our model, although the uniform lithology assumption is used to gauge this deposition-flow-reaction process, different values of the intrinsic permeability in separate gas chimneys and pipe zones have been set to roughly reflect the real lithological situation.One of the two uncertainties lies in the usage of the empirical parameters of SWCC properties in Equations ( 5)-(11).In our simulation, the general values of sedimentary factors are adopted, which define the drainage-imbibition performance of seabed soils.Another is the failure criterion to judge distinct Mohr-Coulomb hydraulic fracturing.We do not elaborate on the scopes of Poisson's ratio, clay content, and the Hoek-Brown constant here, because they are the mean values and can explain most ordinary cases in our research deposits.

Difference of MTDs
In our mode, we do not set any property discrepancies among the three MTDs.The simulation focus is on the lowest MTD3 because it almost marks both the bottom of HSZ and the top of BSR and dominates the unsaturated fluid flow qualities when passing into upper pipe space.Whether or not the overpressure created by upward fluids in gas chimneys gathering beneath MTD3 incites the pipe development via the hydraulic fracturing manner is still controversial.Undoubtedly, the latest emergence period of MTD2 and MTD3 is much earlier than 30 kyr BP because the calculated age of authigenic carbonate rocks at 52.1 mbsf is 0.53 Ma.The authigenic carbonate concretions lying along 3 mbsf level here are estimated to have formed at 3 kyr BP.The shallowest MTD1 is most likely to horizontally insert into the pipe structure.This activity impeded the upward fluid flowwhich just escaped from low MTD2 within the 30 kyr period.That is, multi-stage GH formation might have happened and led to the different occurrences of GHs among MTD1, MTD2, and MTD3 depth layers.

Multiple Accumulation Assumption
The analysis of paleo-GH existence suggests the great possibility that a multi-stage accumulation activity may have occurred.Rapid sedimentation together with sea level variations, landslides, and turbidity currents inevitably induces the simultaneous reaction of GH.The GH formation in the upper and middle layers of the sedimentary layer and decomposition in the lower part occurred [57,59].Huge thick unconsolidated sedimentary seabed breeds and traps large amounts of methane-rich overpressure pore fluids.The mechanical disturbance in early MTD propagation easily triggered the internal rupture of the underlying muddy clay layer and prompted previously isolated pores to connect.The overpressure fluids subsequently penetrated the gas chimney channel, episodically migrated upward, and formed a hydrate-rich layer.After the MTD extension was completed, decreased permeability caused covered layers to serve as a sealing cap.These upward-migrating fluids were impeded.They had to gather largely at the bottom area of the MTDs and develop laterally.Local excellent GH resource sweet-spot layers with a certain thickness and depth were generated in this situation.This effect will be weakened with the temporal and spatial distribution sequence of the MTDs from bottom to top.From the difference in hydrate distribution at site W8, the oldest MTD3 has a GH abundance of up to 50%, and the second oldest MTD2 can also reach 20%.This phenomenon illustrates the rationality of our multiple accumulation speculations.Compared with other highsaturation hydrate reservoirs in oceanic hot reservoirs, such as the Alaminos Canyon in the Gulf of Mexico [60], Ulleung Basin, and the KG Basin offshore Kakinada, their drilling results are often firmly related to the MTD existence.Perhaps it provides a new point of view for hydrate resource exploration optimization.

Conclusions
The GH reservoirs firmly associated with gas chimney, pipe structure, and MTDs have been detected in GMGS5 investigation in the QDN area.How this gas chimney-hydrate geological system evolves and why GH-bearing sediments display the current distribution and occurrence attracts strong interest in probing its dynamic operation mechanism.Because site W08 shows representative system features, it is selected to elucidate the impact of environmental and geological factors, including sea-level change, sedimentation activity, soil-water qualities in and out of the gas chimney area, and MTDs on the evolution process.Through analyzing reasonable deposition-flow-reaction mode matching the local spatiotemporal geological background since 30 kyr BP, six key strata factors at five evolving moments are obtained.The sedimentary tensile fracture and occur at 5 kyr BP.Moreover, the evolution performance of this complex system in six other designed geological settings

Figure 1 .
Figure 1.General geographic information and geological interpretation of the study area: (a) location of QDN area, (b) seafloor slope topography at GMSS5 drilling sites, (c) local seawater bathymetry contours, and (d) comprehensive stratigraphic histogram.

Figure 1 .
Figure 1.General geographic information and geological interpretation of the study area: (a) location of QDN area, (b) seafloor slope topography at GMSS5 drilling sites, (c) local seawater bathymetry contours, and (d) comprehensive stratigraphic histogram.

Figure 3 .
Figure 3. Local depositional conditions including: (a) sea level and sedimentation rate [44], (b) pressure and temperature change at seabed surface, (c) fitting on in situ porosity, and (d) methane concentration from drilling core test and phase lines in different situations (l-liquid, h-GH, and a-gas phase).

Figure 3 .
Figure 3. Local depositional conditions including: (a) sea level and sedimentation rate [44], (b) pressure and temperature change at seabed surface, (c) fitting on in situ porosity, and (d) methane concentration from drilling core test and phase lines in different situations (l-liquid, h-GH, and a-gas phase).

J
. Mar. Sci.Eng.2024, 12, x FOR PEER REVIEW 10 of 19 combination of related key parameter values perfectly explain how gas chimneys and MTD determine the dynamic evolution and final occurrence of GH systems.

Figure 4 .
Figure 4.The evolution state of: (a) the change of pressure difference, (b) the change of temperature difference, (c) gas and hydrate content, (d) pore salinity, (e) capillary pressure, and (f) permeability, at five moments 30, 20, 10, 5 and 0 kyr BP.

Figure 4 .
Figure 4.The evolution state of: (a) the change of pressure difference, (b) the change of temperature difference, (c) gas and hydrate content, (d) pore salinity, (e) capillary pressure, and (f) permeability, at five moments 30, 20, 10, 5 and 0 kyr BP.

Figure 5 .
Figure 5.The evolution state of: (a) tensile strength, and (b) fracturing criterion term Sv at these five moments.

Figure 5 .
Figure 5.The evolution state of: (a) tensile strength, and (b) fracturing criterion term S v at these five moments.

Figure 6 .
Figure 6.The differences of: (a) free gas, (b) GH, and (c) the criterion term Sv between the respective results of these six scenes and the previous basic condition.

Figure 6 .
Figure 6.The differences of: (a) free gas, (b) GH, and (c) the criterion term S v between the respective results of these six scenes and the previous basic condition.

J 19 Figure 7 .
Figure 7.The difference of: (a) free gas, (b) GH, and (c) criterion term Sv with the corresponding value in the basic environment in three cases when initial GH saturation is 5%, 10%, and 15%, respectively.

Figure 7 .
Figure 7.The difference of: (a) free gas, (b) GH, and (c) criterion term S v with the corresponding value in the basic environment in three cases when initial GH saturation is 5%, 10%, and 15%, respectively.

Table 2 .
Six possible geological scenarios for GH generation that may occur in the QDN area.