The study area is located in the deepwater area of the PRMB in the northern SCS ( Figure 1A ). The SCS is the largest marginal sea along the western Pacific and formed as a consequence of multi-phases of continental rifting along the margin of the South China Block from the Late Cretaceous to Paleocene (Holloway, 1982; Briais et al., 1993; Hayes and Nissen, 2005; Ding et al., 2018, 2020). During the Cenozoic era, the SCS underwent various tectonic events, including continental break-up, seafloor spreading, subduction along the Manila Trench, and uplift of the Taiwan Island, accompanying the formation of several sedimentary basins (Yao, 1996; Mi et al., 2018; Zhao et al., 2023a) ( Figure 1A ).

The PRMB lies in the northern continental margin of the SCS, with a total area of 200,000 km² ( Figure 1A ). The basin consists of two depression zones separated by three uplift zones. From north to south, these zones include the northern terrace, the northern depression zone (comprising the Zhu 3 and Zhu 1 depressions), the central uplift zone [comprising the Shenhu Uplift (SHU), Panyu Low Uplift (PYLU), and Dongsha Uplifts (DSU)], the southern depression zone [comprising the Zhu 2 and Chaoshan depressions (CSD)], and the southern uplift zone along the margin of the South China Sea Basin (SCSB) (Xie et al., 2014; Mao et al., 2020) ( Figure 1A ). The formation of the PRMB involved a syn-rift stage during the Eocene epoch (approximately 65 Ma to 32 Ma), followed by a phase of continental breakup in the Oligocene epoch (32 Ma to 23.8 Ma), and concluded with a final post-rift subsidence stage that has persisted since the Miocene epoch (23.8 Ma to the present) (Clift and Lin, 2001; Pang et al., 2008). At 23.8 Ma, a significant tectonic event known as the Baiyun Event occurred, which caused the shelf break zone to shift northward from the southern Baiyun Depression to the PYLU (Pang et al., 2009; Liu et al., 2011). This event also triggered a large-scale transgression in the northern SCS (Pang et al., 2009; Liu et al., 2011). Consequently, the Baiyun sag (BYS), the largest sag in the PRMB, transformed from a shelf environment to its current deep-water slope setting ( Figure 1A ).

The BYS is located in the southern PRMB (Zhou et al., 2009; Zhu et al., 2023) ( Figure 1A ), and appears as NE-trending depression covering an area of over 20,000 km² (Xie et al., 2014). It is surrounded by the PYLU to the north, the Yunli Low Uplift (YLLU) to the south, the SHU to the west, and the DSU to the east (Xie et al., 2014) ( Figure 1A ). Numerous submarine canyons exhibiting characteristics of unidirectional migration have developed along the slope of the northern BYS (Gong et al., 2013; Li et al., 2023) ( Figure 1B ). Furthermore, the BYS is rich in hydrocarbon resources, such as gas hydrate (Yang et al., 2008, 2015, 2017a) ( Figure 1B ). The Wenchang (47.8 to 38 Ma) and Enping formations (38 to 33.9 Ma) serve as the main source rocks characterized by high maturity, which can yield huge amounts of deep thermogenic gas (He et al., 2013; Wang et al., 2018). The Zhujiang (23 to 16 Ma) and Hanjiang formations (16 to 10 Ma) with relatively lower maturity primarily generated biogenic gas (He et al., 2013; Mi, 2018; Wang et al., 2018).

The Shenhu Canyon Group is situated to the north of the BYS (Shenhu area), oriented in NNW-SSE direction ( Figures 1A, B ). The water depths in the Shenhu area range from 500 to 1000 meters, covering an approximate area of 400 km² (Ding et al., 2013).It consists of 19 canyons, labeled as C1 to C19 from west to east along the slope region, which are positioned perpendicular to the continental slope (Li et al., 2023) ( Figure 1B ). The heads of these canyons are confined within the upper continental slope, south of the continental shelf with a width of 236 km (Huang et al., 1995). To the south, the Shenhu Canyon Group intersects with the EW trended Pearl River Canyon, which serves as one of the primary conduits for sediment sourced from onshore areas into deep-water basins (Ding et al., 2013) ( Figure 1B ).

The BYS in the PRMB has been identified as a hydrocarbon-rich basin (Yang et al., 2008; Wang et al., 2011) ( Figure 1A ). Three drilling expeditions have been conducted by the Guangzhou Marine Geologic Survey (GMGS-1, -3, and -4) at the Shenhu Canyon Group ( Figure 1B ). The gas hydrates, investigated during the GMGS-1 drilling expedition in 2007, were sampled at three sites located between C10 and C11 canyons (encompassing sites SH2, SH3, and SH7) (Zhang et al., 2007; Yang et al., 2008) ( Figure 1C ). The GMGS-3 and GMGS-4 expeditions, conducted in 2015 and 2016 respectively, focused in particular in the adjacent region spanning from C12 to C14 canyons. The gas hydrates were sampled at 8 sites (including SC-W01, SC-W02, SC-W03, SH-W07, SH-W11, SH-W17, SH-W18, and SH-W19) (Yang et al., 2015; Yang et al., 2017b) ( Figure 1B ).

The occurrence of BSR is a significant indicator of the presence of gas hydrate (Shedd et al., 2012; Wang et al., 2014). In the Shenhu Canyon Group, previous studies have revealed plenty of BSRs (Chen et al., 2016; Zhang et al., 2017, 2020a; Wang et al., 2021a; He et al., 2022), which are primarily distributed along canyon ridges between canyons C10-C14 and at the toe of canyon C10-C16 ( Figure 1C ). BSRs have been identified within both Pliocene layer and Quaternary layers (He et al., 2022) ( Figure 1C ). However, the presence of BSRs may not guarantee the existence of gas hydrate, as it can also arise from free gas-bearing sediments situated beneath an impermeable layer. Therefore, the identification of BSRs using seismic data must be corroborated with drilling well data to determine the presence of gas hydrate (Yang et al., 2015; Zhang et al., 2020a). Gas chimneys serve as significant pathway for the migration of deep gas to shallower GHSZ, which can provide an adequate gas source for the formation of gas hydrates (He et al., 2022) ( Figure 1C ). BSRs are generally recognized at top of gas chimneys (Zhang et al., 2020b; He et al., 2022). Pockmarks are seafloor notches typically formed by the release of methane gas from subsurface strata through migration channels (Mazzini et al., 2017). The formation of pockmarks is often associated with active fluid overpressure activities, particularly in gas hydrate-rich slope environments (Mazzini et al., 2017). The overpressure fluids that lead to the formation of pockmarks generally originate from the decomposition of gas hydrates (Mazzini et al., 2017). Consequently, pockmarks on the seafloor serve as significant indicators of the decomposition of gas hydrates ( Figure 1C ).

Gas hydrates are stable under conditions of low temperatures and high pressure. To calculate the GHSZ thickness (Miles, 1995), several input parameters are essential. These include water depth, sediment thickness, BWT, and geothermal gradient (Miles, 1995), which has been successfully applied to the neighboring Qiongdongnan Basin (Zhao et al., 2023b) and Dongsha area (Song et al., 2022b). Empirically, a fourth-order polynomial formula is utilized to determine the thickness of GHSZ with a pore-water salinity of 35 wt% and gas source can be assumed to be entire biogenic methane (Miles, 1995),

where a=l.559474×10^-1; b=4.8275×10^-2; c=-2.78083×10^-3 and d=l.5922×10^-4. P_b is the theoretical sediment hydrostatic pressure in MPa, indicating the minimum pressure at which gas and water are combined to form gas hydrates at temperature T in °C.

Sediment temperature is,

where T_z (°C) is the sediment temperature at depth z (meters below sea floor). ΔT / Δz is the geothermal gradient. T₀ is the BWT.

The sediment hydrostatic pressure (MPa) is,

where P_s is the pressure in sediment, D = z₀+z; z is sediment depth below seafloor in m; z₀ is the water depth (m); C₁=(5.92 + 5.25*sin²(Lat))×10^-3; Lat is the latitude (in degrees); C₂ = 2.21×10^-6. Transform Equation 2 into Equation 3, then,

which is the sediment hydrostatic pressure-temperature relationship within the sediment in the same units as in (1). By equating Equations 1, 4, the sediment temperature Tz is calculated. Using Equation 2, we obtain the sediment depth z, which represents the thickness of GHSZ. The sediments are assumed to be isotropic here.

Subsequently, we can calculate the variations of GHSZ thickness during the Quaternary. Sea level fluctuation and BWT variation since 2.5 Ma in the study area are derived from the ODP Site 1148 (Bates et al., 2014) ( Figure 1A ). The sea level fluctuation in the study area since Quaternary is assumed to be equivalent to the change of water depth because of only 0-10 m tectonic subsidence occurred since Quaternary in the study area (Zhao et al., 2011; Xie et al., 2014). We also assume that the whole study area has the same BWT variation as the ODP site 1148 because our study area is relatively adjacent to this site. The sediment thickness variations since 2.5 Ma are reconstructed based on a sedimentation rate of 9.8 cm/ka at the LW3-1 well from the Pliocene to the Quaternary (Sun et al., 2022; Xie et al., 2021). Based on these data, the thickness variations of GHSZ in the whole study area since 2.5 Ma were calculated ( Figure 3 ).

The thickness of the present GHSZ calculated by using water depth, BWT, geothermal gradient and sediment thickness parameters is approximately 0 m at a water depth of 600 m, to the north of the BYS. Subsequently, it gradually increases, reaching a maximum thickness of 350 m in the YLLU. After that, the GHSZ thickness thins once more in the Liwan Sag (LWS) ( Figure 3 ). Generally, the GHSZ of > 300-m-thick is mainly situated in the south of Shenhu, the southern BYS, the YLLU and the northern LWS ( Figure 3 ).

The drilling well SH-W18 is located in the Shenhu area at a water depth of 1288 m below sea level (mbsl) where BWT is 4.84°C ( Figure 3 ). The thickness of present GHSZ calculated using the Miles (1995) method is 161 m at drilling well SH-W18. On the basis of ODP site 1148, the reconstructed sea level and BWT variation since 2.5 Ma at well SH-W18 ( Figure 4 ) shows that sea level rises or falls are synchronously with the changes of BWT ( Figure 4A ). However, the variation of sea level and BWT have opposite effects on the GHSZ thickness. As BWT decreases, GHSZ thickness increases, whereas as sea level drops, GHSZ thickness decreases. Due to differences in the magnitude of sea level and BWT changes since 2.5 Ma ( Figure 4A ), reduction in the GHSZ thickness is mainly controlled by three patterns of the variation of the sea level fluctuations and BWT ( Figure 4 ). During the stage from 2.5 to 0.9 Ma, the reduction in GHSZ thickness was characterized by a slight rise in sea level and a rapid increase in the BWT ( Figure 4A ). The variation trend of GHSZ thickness caused by BWT change alone with the sea level staying constant at 2.5Ma extremely consistent with that caused by the combined changes of BWT and sea level, and the error is generally within 5m ( Figure 4C ). However, the variation trend of GHSZ thickness caused by sea level change alone with the BWT staying constant at 2.5Ma is significantly different from that caused by the combined changes of BWT and sea level ( Figure 4C ). During the stage from 0.9 to 0.02 Ma, the reduction in GHSZ thickness was characterized by a rapid drop in sea level and a slight decline in the BWT ( Figure 4A ). Compared with the variation trend of GHSZ thickness caused by BWT change alone with the sea level staying constant at 2.5Ma, the variation trend of GHSZ thickness caused by sea level change alone with the BWT staying constant at 2.5Ma is more consistent with the variation trend caused by the combined changes of BWT and sea level, and the error is generally within 5m ( Figure 4C ). During the stage from 0.02 to 0 Ma, the reduction in GHSZ thickness was characterized by a rapid increase both in the sea level and BWT ( Figure 4A ). The variation trend of the GHSZ thickness caused either by the sole change in BWT or by the sole change in sea level is significantly different from that caused by the combined changes of both ( Figure 4C ). Consequently, we selected three periods 2.27-2.26 Ma (increase by 3 m in sea level and 2°C in the BWT), 0.31-0.14 Ma (decrease by 100 m in sea level and 0.5°C in the BWT) and 0.02-0 Ma (increase by 120 m in sea level and 2.6°C in the BWT) (referred to as patterns 1 to 3 here) representing the most significant changes among the three stages to analyze how the GHSZ reduction distribute across the whole study area ( Figures 5 – 7 ).

The BGHSZ is generally regarded as the bottom of the GHSZ, where gas hydrates can exist stably. The BGHSZ is generally visible in seismic data as a BSR (Yang et al., 2021). Based on the water depth, in situ BWT and geothermal gradient data from wells Sh2, Sh3, SC-W01 and SH-W18 (Su et al., 2016; Song et al., 2022a), the GHSZ thickness is calculated by using Miles (1995) method. Subsequently, the BGHSZ was determined and delineated on the seismic profile through the application of a time-depth conversion equation (Su et al., 2018; Wang et al., 2021a) ( Figure 8 ). At Site Sh2, the BSR is located at 221 mbsf (meter below sea floor) based on LWD (Su et al., 2016) with the BGHSZ at 215 mbsf ( Figure 8B ). For Site Sh3, the BSR is at 196 mbsf (Su et al., 2016) with the BGHSZ at 190 mbsf ( Figure 8B ). At SC-W01, the BSR is at 175 mbsf (Song et al., 2022a) with the BGHSZ at 154 mbsf ( Figure 8C ). Similarly, at SC-W18, the BSR is at 172 mbsf (Zhang et al., 2020a) while the modeled BGHSZ (assuming 100% CH4) at well SH-W18 is located at 161 mbsf ( Figure 8C ), which is consistent with the results calculated by Zhang et al. (2020a) using the Sloan method (Sloan, 1998). By comparing of the depths of BSR and BGHSZ at the drilling wells, the results demonstrate that the calculated depths of the BGHSZ are strikingly similar (error generally less than 7%) to the seismically interpreted BSRs ( Figure 8D ). However, the depths of the seismically interpreted BSRs are constantly slightly deeper than that of the modeled BGHSZ.

The gas composition of the study area is important to the numerical calculation for GHSZ thickness. Geochemical analyses of pore water obtained from core samples by GMGS1 expedition suggest that methane concentrations ranged from 98.90% to 99.96%, with minor quantities of ethane and propane (Wu et al., 2008; Su et al., 2016). These scholars believe that gas hydrates are mainly derived from biogenic methane (Wu et al., 2011; Zhu et al., 2013). Geochemical analysis of core samples by GMGS-3 and 4 indicates that the gas composition was predominantly methane gas. Heavy hydrocarbons, including ethane, propane, were also detected in some samples (Wei et al., 2018; Zhang et al., 2019); however, the content of heavy hydrocarbons is at a relatively low level. These scholars contend that gas hydrates are mainly derived from thermogenic methane (Zhang et al., 2019). The deeper depths of seismically interpreted BSRs than the modeled BGHSZ may be explained by the presence of heavy hydrocarbons in the gas composition at these stations, which may potentially cause the BSR to shift downward (Zhang et al., 2020a). Zhang et al. (2019) showed that the modeled BGHSZ (assuming 100% CH₄ as Structure I gas hydrates) is shallower than BSR, the depth of which is consistent with the modeled BGHSZ (assuming 99% CH₄ and 1% C₂H₆ as Structure II gas hydrate). The Result suggested both Structure I and II gas hydrate may coexist in the Shenhu area (Zhang et al., 2019). Wei et al. (2018) confirmed coexistence of the Structure I and II gas hydrate at well SC-W01 in the Shenhu area using Raman spectra analysis based on LWD data and gas hydrate samples obtained by GMGS-4 and suggested that existence of Structure II gas hydrate potentially deepen the base of Structure I gas hydrate (Wei et al., 2018). Qian et al. (2018) reported that Structure II gas hydrate and free gas are interbedded and can coexists below Structure I gas hydrate using well log and core data at well SH-W07 (Qian et al., 2018) ( Figure 1C ). Consequently, our results of the GHSZ are reasonably sound, with the error generally within 7%. The primary cause of the error is the heterogeneity of natural gas components in the study area.

As far as we know, the GHSZ is primarily controlled by temperature and pressure conditions (Chen et al., 2019; Ketzer et al., 2020; Davies et al., 2021; Deng et al., 2021; Li et al., 2024). However, our results show that the reduction in the thickness of the GHSZ in shallow water areas is always more significant than that in deep water areas ( Figures 5 - 7 ). Thus, a series of tests are meticulously carried out to investigate the variations in the thickness of the GHSZ in relation to the alterations in sea level and BWT at diverse water depths ( Figures 11 ).

In Figure 11A , the higher the BWT, the greater the reduction in the thickness of the GHSZ. Moreover, in shallow water areas, especially those with a water depth of less than 1000m, the reduction in the thickness of the GHSZ is significantly more pronounced than that in deep water areas. In Figure 11B , for the same amount of water depth change (for example, an increase of 100m), the change in the thickness of the GHSZ in shallow water areas is always greater than that in deep water areas. The thickness of the GHSZ varies with changes in BWT and sea level. However, in shallow water regions, particularly those with a water depth less than 1000m, the variation in the thickness of the GHSZ is consistently more evident than that in deep water regions. As indicated by the slopes of the lines in Figure 11 , The slopes of the lines in shallow water areas, especially those with a water depth of less than 1000m, are always greater than that in deep water areas.

The Shenhu area is characterized by a shallow water depth ranging from 500 to 1000m and a relatively high current BWT of 3-5 °C. In contrast, the YLLU and LWS areas feature a deep water depth exceeding 1500m and a relatively low current BWT below 3°C ( Figures 2A, B ). The shallow water depths (< 1000m) always have a larger effect on the GHSZ thickness variation than the deep water depths (>1000m) and that is why the Shenhu area always shows a larger reduction in the GHSZ thickness compared to the YLLU and LWS ( Figures 5 - 7 ).

Reduction in GHSZ thickness may lead to gas hydrate decomposition accompanied by the seepage of methane (Marín-Moreno et al., 2013; Ketzer et al., 2020; Davies et al., 2021; Deng et al., 2021; Cooke et al., 2023; Li et al., 2024). Authigenic carbonates were revealed on the seafloor in the Shenhu area (Qiu et al., 2018; Chen et al., 2022; Zhang et al., 2022). The enriched ¹⁸O suggests that the presence of authigenic carbonate is linked to the decomposition of gas hydrates, a process accompanied by methane seepage (Tong et al., 2013). This phenomenon is also observed in the Dongsha area and the Qiongdongnan Basin (Zhang et al., 2023; Li et al., 2024) ( Figure 4A ). Methane seepage can breach the seafloor and form seafloor pockmarks (Mazzini et al., 2017; Qiu et al., 2018). Observations of seismic profiles reveal seafloor pockmarks are presented on both sides of the Shenhu Canyon Group, with some situated above the BSR, suggesting that gas hydrates may have decomposed or be currently decomposing (Chen et al., 2016) ( Figures 9B , 10C ). Paleo-pockmarks have also been observed in seismic profiles ( Figure 10C ), suggesting that fluid seepage activity occurred in the past (Wang et al., 2021b; Shen et al., 2023). So far, the history of methane leakage in the northern SCS has been recorded for only 0.4 Ma (Tong et al., 2013; Qiu et al., 2018; Chen et al., 2019; Deng et al., 2021; Chen et al., 2022; Zhang et al., 2022, 2023; Li et al., 2024) ( Figure 4A ). During this period, methane leakage occurred either during sea level falls ( Figure 6 ) or BWT increase ( Figures 5 , 7 ), which is consistent with our findings ( Figures 4 - 7 ). Consequently, methane leakage occured in the Shenhu area due to the decomposition of gas hydrates. Apart from the influence of the BWT and sea level fluctuations, water depth is a crucial factor affecting gas hydrate decomposition. In shallow water areas, especially those with a water depth of less than 1000 m, gas hydrate decomposition is consistently more significant than that in deep water areas.

Our findings reveal significant fluctuations in the GHSZ thickness since 2.5 Ma, with multiple thinning periods controlled by the above-mentioned three combinations of sea level and BWT variations ( Figures 5 - 7 ). Due to predominant influence of shallow water area on the variation of GHSZ thickness ( Figure 11 ), the reduction of GHSZ thickness, contributing to gas hydrates leakage, have occurred frequently in the Shenhu area ( Figures 4 ). Evidences of double BSRs, free gas, seafloor pockmarks, paleo-pockmark in Seismic data ( Figures 10 – 12 ) and authigenic carbonates (Tong et al., 2013; Chen et al., 2022; Zhang et al., 2022) prove that gas hydrate leakage has occurred in the Shenhu area.

Gas hydrate dissociation is known to trigger methane leakage and form pockmarks (Chen et al., 2019; Ketzer et al., 2020; Davies et al., 2021; Deng et al., 2021; Chen et al., 2022; Zhang et al., 2022; Cooke et al., 2023; Panieri et al., 2023; Zhang et al., 2023; Li et al., 2024). The leakage can increase pore pressure within sediment layers and destabilize the slope, ultimately resulting in slope failures (Sultan et al., 2004; Mienert et al., 2005; Nixon and Grozic, 2006, 2007). Slope failures were observed in both sides of canyon walls ( Figure 9B ). Yu et al. (2023) carried out numerical simulations to reconstruct the erosional process of pockmarks induced by bottom currents in the northwestern SCS. The simulation results clearly demonstrated that the decomposition of gas hydrates and the seepage of methane can result in the formation of pockmarks on the seafloor. Moreover, specific bottom current velocities are capable of eroding these pockmarks, causing them to connect with one another and form trains of pockmarks. Eventually, a canyon will be formed.

Thus, we propose a possible formation mechanism for the Shenhu Canyon Group, which is as follows: