Geofluids
Volume 2019 (2019), Article ID 1010824, 18 pages
https://doi.org/10.1155/2019/1010824

Research Article

Methane Source and Turnover in the Shallow Sediments to the West of Haima Cold Seeps on the Northwestern Slope of the South China Sea

Junxi Feng,1,2 Shengxiong Yang,¹ Hongbin Wang,¹ Jinqiang Liang,¹ Yunxin Fang,¹ and Min Luo3,4

¹MLR Key Laboratory of Marine Mineral Resources, Guangzhou Marine Geological Survey, Guangzhou 510075, China
²School of Marine Sciences, Sun Yat-sen University, Guangzhou 510006, China
³Shanghai Engineering Research Center of Hadal Science and Technology, College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China
⁴Laboratory for Marine Geology, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266061, China

Correspondence should be addressed to Shengxiong Yang; yangshengxiong@gmgs.cn and Min Luo; mluo@shou.edu.cn

Received 16 February 2019; Revised 7 May 2019; Accepted 27 May 2019; Published 5 August 2019

Academic Editor: Giovanni Martinelli

Copyright © 2019 Junxi Feng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The Haima cold seeps are active cold seep areas that were recently discovered on the northwestern slope of the South China Sea (SCS). Three piston cores (CL30, CL44, and CL47) were collected within an area characterized by bottom simulating reflectors to the west of Haima cold seeps. Porewater profiles of the three cores exhibit typical kink-type feature, which is attributed to elevated methane flux (CL30) and bubble irrigation (CL44 and CL47). By simulating the porewater profiles of SO₄²⁻, CH₄, PO₄³⁻, Ca²⁺, Mg²⁺, and dissolved inorganic carbon (DIC) in CL44 and CL47 using a steady-state reaction-transport model, we estimated that the dissolved SO₄²⁻ was predominantly consumed by anaerobic oxidation of methane (AOM) at rates of 74.3 mmol m⁻² yr⁻¹ in CL44 and 85.0 mmol m⁻² yr⁻¹ in CL47. The relatively high AOM rates were sustained by free gas dissolution rather than local methanogenesis. Based on the diffusive Ba²⁺ fluxes and the excess barium contents in the sediments slightly above the current SMTZ, we estimated that methane fluxes at core CL44 and CL47 have persisted for ca. 3 kyr and 0.8-1.6 kyr, respectively. The non-steady-state modeling for CL30 predicted that a recent increase in upward dissolved methane flux was initiated ca. 85 yr ago. However, the required time for the formation of the barium front above the SMTZ at this core is much longer (ca. 2.2-4.2 kyr), which suggests that the depth of SMTZ possibly has fluctuated due to episodic changes in methane flux. Furthermore, using the model-derived fractions of different DIC sources and the δ¹³CDIC mass balance calculation, we estimated that the δ¹³C values of the external methane in cores CL30, CL44, and CL47 are -74.1‰, -75.4‰, and -66.7‰, respectively, indicating the microbial origin of methane. Our results suggest that methane seepage in the broader area surrounding the Haima cold seeps probably has persisted at least hundreds to thousands of years with changing methane fluxes.

1. Introduction

Methane in marine sediments as dissolved gas in porewater or free gas (bubbles) depending on its in situ solubility is a significant component of the global carbon cycle. Methane could also exist in an ice-like solid as gas hydrate if the in situ gas hydrate solubility concentration is oversaturated at suitable pressure-temperature conditions [1]. The base of the gas hydrate reservoir in marine sediments is present as a characteristic discontinuity known as a bottom-simulating reflector (BSR), which results from the occurrence of free gas beneath the gas hydrate stability zone (GHSZ) [2]. The great majority of methane is consumed by the microbial consortium via anaerobic oxidation of methane within the sulfate-methane transition zone (SMTZ) where methane meets sulfate diffusing downwards from seawater (AOM: ) [3, 4]. Through this, reaction methane is converted to dissolved inorganic carbon (DIC) which could be partially removed from solution by authigenic carbonate precipitation [5, 6]. Therefore, AOM largely prevents dissolved methane from entering water column and plays a significant role in marine carbon cycling.

Gas bubble rise is a particularly effective mechanism for transporting methane through the sediment and into the bottom water because gas ascension can be faster than bubble dissolution [7] and methane gas cannot directly be consumed by microorganisms [8]. The rising methane gas bubbles emitting from the seafloor can mix bottom seawater down into the sediment column over several meters. The resulting kink-type porewater profiles are supposed to be stable for several years to decades even after active gas bubble ebullition has ceased [7]. Enhancement in upward methane flux could also result in kink-type or concave-up porewater profiles [9–12]. Hence, these types of nonlinear porewater profiles can thus be used to estimate the timing of (sub)recent methane pulse and provide insights into the dynamics of methane seepage and the underlying gas hydrate reservoir [10, 11, 13, 14].

At the steady-state condition, Ba²⁺ diffusing upward into the sulfate-bearing zone above the SMTZ precipitates as barite and forms the authigenic barium fronts () [15]. When buried beneath the sulfate-bearing zone, barite tends to dissolve and release Ba²⁺ into the porewater below the SMTZ due to unsaturation. Through this cycling, authigenic barium fronts would be stably developed just above the SMTZ [15–17]. The content of authigenic barite depends on the upward diffusive Ba²⁺ flux and the duration that the SMTZ has persisted at a given depth interval. The time required for barium front formation above the SMTZ could thus be calculated based on the depth-integrated excess barium contents and the porewater dissolved barium concentration gradients assuming a constant upward Ba²⁺ flux. Therefore, the authigenic barium fronts in sediments can be used to trace present and past SMTZ and associated methane release events as well as the duration of methane seepage that has persisted under a given methane flux [16–20].

Methane seepages are widespread on the northern slope of the South China Sea (SCS) as revealed by authigenic carbonates collected at more than 30 cold seep sites [21–26]. The Haima cold seeps were recently discovered on the northwestern slope of SCS [25]. Several sites with gas bubbling identified by hydroacoustic anomalies, and shallow gas hydrates were found around this area [27–31]. Recent studies have shown a pronounced temporal change in methane seepages and a potential lateral migration of methane-bearing fluid along more permeable sand-bearing layer at Haima cold seeps [25, 32, 33]. Nevertheless, our quantitative understanding of the methane dynamics in this area remains scarce.

In this study, we present porewater geochemical data of three piston cores (CL30, CL44, and CL47) collected to the west of Haima cold seeps, including concentrations of sulfate (SO₄²⁻), calcium (Ca²⁺), magnesium (Mg²⁺), barium (Ba²⁺), phosphate (PO₄³⁻), methane (CH₄), and DIC as well as the carbon isotopic compositions of DIC. Using a steady-state reaction-transport model, we quantify the methane turnover rates in CL44 and CL47 which are mainly supplied by rising free gas. The kink in the porewater profiles of CL30 was reconstructed using a non-steady-state modeling approach assuming a recent increase in methane flux. In addition, authigenic Ba enrichments were used to constrain the durations that the current or past methane seepages have persisted. Furthermore, a simple mass balance model of DIC and δ¹³CDIC was applied to explore the methane source.

2. Geological Background

The northern SCS is characterized as a Cenozoic, Atlantic-type passive continental margin [34], where the marginal basins generally underwent two stages of evolution, including the rift stage and the postrift thermal subsidence stage [35]. Qiongdongnan Basin is a northeastern trended Cenozoic sedimentary basin which developed on the northwestern part of the SCS [36]. Covered by sedimentary materials of up to 10 km, the depositional environment of the basin initially transformed from lacustrine to marine conditions and later from neritic to bathyal, starting from Eocene till present [37]. During the rifting stage, numerous half-grabens and sags were developed. After that, postrift thermal subsidence occurred and a thick sediment sequence dominated by mudstones was deposited in the basin since Miocene. Collectively, the sedimentation rates and the present-day geothermal gradient are both high in the Qiongdongnan Basin [38]. The thick sediment sequences, high geothermal gradient along with faulting and/or diapirism, have facilitated the generation and migration of the hydrocarbons in the basin [39]. The widely distributed bottom-simulating reflectors and gas chimneys identified in the Qiongdongnan Basin were linked to the accumulation of gas hydrate [40, 41].

The active Haima cold seeps have been discovered in the southern uplift belt of the Qiongdongnan Basin on the lower continental slope of the northwestern SCS during R/V Haiyang-6 cruises in 2015 and 2016. Abundant chemosynthetic communities, methane-derived authigenic carbonates, and massive gas hydrates were found at the Haima cold seeps [25]. The dating of bivalve shells and seep carbonates revealed episodic changes in seepage activity [25]. Other features of methane seeps, such as acoustic plume, acoustic void, chimney structures, and pockmarks, were also reported at the Haima cold seeps and its surrounding area [27–31, 42]. The sampling sites are ca. 20 to 30 kilometers west of the Haima cold seeps, where BSR is well developed (Fang Y., unpublished data). The bathymetric investigation has shown a relatively flat topography, and the water depths range from 1250 to 1300 m in the study area (Figure 1).

(a)

(b)

Figure 1: (a) Location of the study area. The grey area represents the subsurface area lineated by seismic investigation. The location of the reference core SO49-37KL (blue circle) is also shown. (b) Locations of sampling sites (red dots) and the Haima cold seeps (blue dots).

3. Materials and Methods

3.1. Sampling and Analytical Methods

Three piston cores (CL30, CL44, and CL47) were collected from the southern Qiongdongnan Basin west to the Haima cold seeps at water depths ranging from 1255 m to 1301 m during the R/V Haiyang-4 cruise conducted by Guangzhou Marine Geological Survey in 2014 (Figure 1 and Table 1). The sediments of the three cores mainly consist of greyish-green silty clay. Notably, the sediments at the bottom of core CL44 yielded a strong odour of hydrogen sulfide. Porewater samples were then collected onboard using Rhizon samplers with pore sizes of the porous part of approximately 0.2 mm at intervals of 20 cm for CL44 and 60 cm for CL30 and CL47. All the porewater samples were preserved at ~4°C until further analyses.

Table 1: Information on the studied cores from the northwestern South China Sea.

Site
Water depth (m)
Seafloor temperature (°C)
Core length (cm)

CL30
1255
3.3
630

CL44
1279
3.2
752

CL47
1301
3.1
775

PO₄³⁻ concentrations were measured onboard using the spectrophotometric method according to Grasshoff et al. [43] with a UV-Vis spectrophotometer (Hitachi U5100). The precision for phosphate was ±3.0%. 10 ml of sediments was added to 20 ml empty vials onboard to replace the 10 ml headspace needed for the chromatograph injection. The concentrations of hydrocarbon gas were measured onboard using the gas chromatograph method (Agilent 7890N). The precision for methane measurements was ±2.5% [44]. Porosity and density were determined from the weight loss before and after freeze-drying of the wet sediments using a cutting ring with definite mass (15 g) and volume (9.82 cm³) onboard at core CL44. The porosity and density were calculated assuming a density of the porewater of 1.0 g cm⁻³.

The offshore analyses of porewater samples for core CL44 and for cores CL30 and CL47 were performed at the Nanjing University and the Third Institute of Oceanography, State Oceanic Administration, respectively. For core CL44, SO₄²⁻, Ca²⁺, and Mg²⁺ were measured using the standard method of ion chromatography (Metrohm 790-1, Metrosep A Supp 4-250/Metrosep C 2-150). The relative standard deviation was less than 3%. Ba²⁺ concentrations were measured by inductively coupled plasma mass spectrometry (ICP-MS, Finnigan Element II). Before measurement, samples were prepared by diluting in 2% HNO₃ with 10 ppb of Rh as an internal standard. The analytical precisions were estimated to be <5% for Ba²⁺. For cores CL30 and CL47, SO₄²⁻, Ca²⁺, and Mg²⁺ concentrations were determined on a Thermo Dionex ICS-1100 ion chromatograph after a 500-fold dilution using ultrapure water [44]. Porewater samples were prepared by diluting in 2% HNO₃ with 10 ppb of Tb as an internal standard before analysis for Ba²⁺ using the ICP-MS (Thermo Fisher iCAPQ). The analytical precisions were estimated to be <5% for Ba²⁺.

For core CL44, DIC concentrations and δ¹³CDIC values were determined using a continuous flow mass spectrometer (Thermo Fisher Delta-Plus). 0.5 ml porewater was treated with pure H₃PO₄ in a glass vial at 25°C. The CO₂ produced was stripped with He and transferred into the mass spectrometer through the measurement of the δ¹³C value [45]. For cores CL30 and CL47, the DIC concentrations and carbon isotopic ratios were determined via a continuous flow mass spectrometer (Thermo Delta V Advantage). A 0.2 ml porewater sample was treated with pure H₃PO₄ in a glass vial at 25°C. The CO₂ produced was stripped with He and transferred into the mass spectrometer through which the δ¹³C values were measured. The analytical precisions were better than 0.2‰ for δ¹³C and better than 2% for DIC concentration [44].

The particulate organic carbon (POC) contents were determined using the potassium dichromate wet oxidation method. The relative standard deviation of the POC content is <1.5%. The aluminium (Al), silicon (Si), and titanium (Ti) concentrations of the sediment samples at cores CL30, CL44, and CL47 were analyzed using PANalytical AXIOSX X-ray fluorescence spectrometry (XRF). The analytical precisions were estimated to be <2% for Al, Si, and Ti. The contents of Ba, zircon (Zr), and rubidium (Rb) in bulk sediments were determined using a PerkinElmer Optima 4300DV ICP-OES after digestion using HCl, HF, and HClO₄ acid mixture. Rhodium was added as an internal standard for calculating the concentrations of the trace elements. The analytical precisions were estimated to be <2% for Ba, Zr, and Rb. The carbonate (CaCO₃) contents of the sediment samples were determined by titration with EDTA standard solution. The analytical precisions were estimated to be <2%. For grain size measurements, approximately 0.5 g of the unground sample was treated with 10% () H₂O₂ for 48 h to oxidize organic matter and then dispersed and homogenized in sodium hexametaphosphate solution using ultrasonic vibration for 30 s before being analyzed by a laser grain size analyzer (Mastersizer 2000). The detection limit ranged from 0.5 μm to 2000 μm. in size were classified as clay, 4 to 63 μm as silt, and larger than 63 μm as sand. The analytical precision is better than 3%.

3.2. Diffusive Flux Calculation

To calculate the diffusive Ba²⁺ fluxes below the kink at cores CL30, CL44, and CL47, equations (1) and (2) were used assuming a steady-state condition [46]:

where represents the diffusive flux of Ba²⁺ (mmol m⁻² yr⁻¹), is the porosity, is the diffusion coefficient for seawater (m² s⁻¹), is the diffusion coefficient for sediments (m² s⁻¹), is the concentration of barium (mmol l⁻¹), and x is the sediment depth (m). The average of sediment porosity of core CL44 (0.69) is applied to cores CL30 and CL47.

3.3. Estimating the Accumulation Time of Diagenetic Barite

The total amount of excess Ba within the interval of the barium peak was calculated using an integral equation:

where is the integral value of barium concentration in a peak from a depth interval from to , and are the average grain density and porosity of the sediments, respectively.

Under the premise of a constant diffusive upward flux of Ba²⁺ into the sulfate-bearing zone, the time needed for barium front formation was calculated using the equation:

In this case, is the time for barite enrichment, is the depth-integrated excess barium content within a peak, and is the upward diffusive flux of Ba²⁺. The diffusive flux was calculated using equations (1) and (2). is the tortuosity- and temperature-corrected diffusion coefficient of Ba²⁺ in the sediment, calculated from the diffusion coefficient in free solution () of 4.64, 4.62, and  cm² s⁻¹ (3.3, 3.2, and 3.1°C) for CL30, CL44, and CL47, respectively, according to Boudreau [47].

3.4. Reaction-Transport Model

A one-dimensional, steady-state, and reaction-transport model was applied to simulate one solid (POC) and six dissolved species including SO₄²⁻, CH₄, DIC, PO₄³⁻, Ca²⁺, and Mg²⁺. The model is modified from previous simulations of methane-rich sediments [48–51], and a full description of the model is shown in Supplementary Materials. All the reactions considered in the model and the expression of kinetic rate are listed in Table 2.

Table 2: Rate expressions of the reactions considered in the model.

Rate
Kinetic rate law

Total POC degradation (wt.% C yr⁻¹)

POM degradation via sulfate reduction (mmol cm⁻³ yr⁻¹ of SO₄²⁻)

Methanogenesis (mmol cm⁻³ yr⁻¹ of CH₄)

Anaerobic oxidation of methane (mmol cm⁻³ yr⁻¹ of CH₄)

Authigenic Ca-carbonate precipitation (mmol cm⁻³ yr⁻¹ of Ca²⁺)

Authigenic Mg-carbonate precipitation (mmol cm⁻³ yr⁻¹ of Mg²⁺)

Gas bubble irrigation (mmol cm⁻³ yr⁻¹)

Gas bubble dissolution (mmol cm⁻³ yr⁻¹ of CH₄)

fPOC converts between POC (dry wt.%) and DIC (mmol cm⁻³ of porewater): , where is the molecular weight of carbon (12 g mol⁻¹), is the density of dry sediments, and is the porosity.

Solid species are transported through the sediments only by burial with prescribed compaction, which is justified because we are only concerned with the anoxic diagenesis below the bioturbated zone. For sites CL44 and CL47, solutes are considered to be transported by molecular diffusion, porewater burial, and gas bubble irrigation, whereas for site CL30, solutes are regarded to be transported by molecular diffusion and porewater burial. Rising gas bubbles facilitate the exchange of porewater and bottom water as they move through tube structures in soft sediments [7]. Although this process was not observed directly, there are evidences implying that it is a significant pathway for transporting methane into the upper 10 m of sediment at sites CL44 and CL47 and driving the mixture of porewater and seawater in the upper two meters (see Section 5.1). The induced porewater mixing process was described as a nonlocal transport mechanism whose rate for each species is proportional to the difference between solute concentrations at the sediment surface (mmol cm⁻³) and at depth below the sediment surface (mmol cm⁻³) (, Table 2). Bubble irrigation is described by parameters (yr⁻¹) and (cm) that define the irrigation intensity and its attenuation below the irrigation depth (cm), respectively [49]. The latter can be determined by visual inspection of the porewater data (see Results) whereas is a model fitting parameter. For the sake of parsimony, is assumed to be constant for both sites.

Although dissolution of gas was allowed to occur over the whole sediment column, the rising methane gas was not explicitly modeled. The rate of gas dissolution, (mmol cm⁻³ yr⁻¹), was described using a pseudo-first-order kinetic expression of the departure from the local methane gas solubility concentration, (mmol cm⁻³), where (yr⁻¹) is the kinetic constant for gas bubble dissolution (Table 2). Methane only dissolves if the porewater is undersaturated with respect to :

was calculated for the in situ salinity, temperature, and pressure using the algorithm in [52]. was constrained using the dissolved sulfate and DIC data (see below).

Major biogeochemical reactions considered in the model are particulate organic matter (POM) degradation via sulfate reduction, methanogenesis, AOM, and authigenic carbonate precipitation. Organic matter mineralization via aerobic respiration, denitrification, and metal oxide reduction were ignored since these processes mainly occur in the surface sediments which were mostly lost during coring.

POM is chemically defined as , where CH₂O and POP denote particulate organic carbon and phosphate, respectively. The total rate of POM mineralization, (wt.% C yr⁻¹), is calculated by the power law model from [53] that considers the initial age of organic matter in surface sediments, (yr) (Table S2). POM mineralization coupled to sulfate reduction follows the stoichiometry:

where is the ratios of particulate organic phosphate to carbon. It is assumed to be the typical ratios as 1/106 [48].

When sulfate is almost completely consumed, the remaining POM is degraded via methanogenesis:

The dominant pathways of methanogenesis in marine sediments are organic matter fermentation and CO₂ reduction [54]. Their net reactions at steady state are balanced with equivalent amounts of CO₂ and CH₄ being produced per mole of POM degraded [55]. Therefore, the reaction of methanogenesis is a net reaction.

Methane is considered to be consumed by AOM [3]:

The rate constant for AOM, (cm³ mmol⁻¹ yr⁻¹), is tuned to the sulfate profiles within the SMTZ.

The loss of Ca²⁺ and Mg²⁺ resulting from the precipitation of authigenic carbonates as Ca-calcite and Mg-calcite was simulated in the model using the thermodynamic solubility constant as defined in [56] (Table 2). A typical porewater pH value of 7.6 was used to calculate CO₃²⁻ from modeled DIC concentrations [57]. was not simulated explicitly in the model.

The length of the simulated model domain was set to 1000 cm. Upper boundary conditions for all species were imposed as fixed concentrations (Dirichlet boundary) using measured values in the uppermost sediment layer where available. For CL44 and CL47, a zero concentration gradient (Neumann-type boundary) was imposed at the lower boundary for all the species. For CL30, a zero concentration gradient was imposed at the lower boundary for all the species except CH₄. CH₄ concentration at the lower boundary was a tunable parameter constrained from the SO₄²⁻ profile. The model was solved using the NDSolve object of MATHEMATICA V. 10.0. The steady-state simulations were run for 10⁷ yrs to achieve the steady state with a mass conservation of >99%. Further details on the model solutions can be found in Supplementary Materials. For the non-steady-state modeling of CL30, a fixed methane concentration in equilibrium with the gas hydrate solubility constrained by local seafloor temperature, pressure, and salinity was defined as the lower boundary of methane [58]. The extrapolation of sulfate concentrations in the upper 3.5 m to zero was taken as the initial condition prior to the increase in methane flux (Supplementary Materials). The basic model construction and kinetic rate expressions as well as the upper and lower boundary conditions for other species were identical to those in the steady-state model.

4. Results

4.1. General Geochemical Trends

The depth profiles of SO₄²⁻ concentration showed kink-type features at all the three cores (Figures 2–4 and Table 3). At site CL30, SO₄²⁻ concentrations decreased gradually above a kink at ~3.5 mbsf and the gradient became steeper below that depth towards the SMTZ at ~4.7 mbsf (Figure 2). In contrast, SO₄²⁻ concentrations at sites CL44 and CL47 displayed near-seawater values in the upper ~2 mbsf above the kinks and then decreased sharply down to the SMTZ located at ~7 and ~6.8 mbsf, respectively (Figures 3 and 4). Ca²⁺ and Mg²⁺ concentrations showed similar trends, with gradual decrease in the upper layers at core CL30 and close to seawater concentration above the kinks at cores CL44 and CL47. Ca²⁺ and Mg²⁺ concentrations declined sharply below the kinks due to ongoing carbonate precipitation and reached minimum at the SMTZ (Figures 2–4 and Table 3). Concentrations of DIC and PO₄³⁻ showed opposite trends to SO₄²⁻, being depleted within the upper layer and enriched below it with the maximum at the STMZ (Figures 2–4 and Table 3). Moreover, CH₄ concentrations at the three cores sharply increased below the SMTZ. The scatter in the CH₄ contents was due to the degassing during core retrieval. The DIC concentrations increased with depth and reached maximum at the SMTZ, with opposite trends of δ¹³CDIC values (minimum values: -46.4‰ for CL30, -41.0‰ for CL44, and -38.8‰ for CL47) (Figures 2–4 and Table 3).

Figure 2: Measured (dots) and simulated (curves) depth profiles of core CL30. Down-depth concentration of particulate organic carbon (POC), sulfate (SO₄²⁻), methane (CH₄), phosphate (PO₄³⁻), dissolved inorganic carbon (DIC), calcium (Ca²⁺), magnesium (Mg²⁺), and δ¹³CDIC is shown.

Figure 3: Measured (dots) and simulated (curves) depth profiles of core CL44. Down-depth concentration of POC, SO₄²⁻, CH₄, PO₄³⁻, DIC, Ca²⁺, Mg²⁺, and δ¹³CDIC is shown.

Figure 4: Measured (dots) and simulated (curves) depth profiles of core CL47. Down-depth concentration of POC, SO₄²⁻, CH₄, PO₄³⁻, DIC, Ca²⁺, Mg²⁺, and δ¹³CDIC is shown.

Table 3: Concentrations and isotope ratios of various dissolved components at cores CL30, CL44, and CL47.

Depth (cmbsf)
CH₄ (mM)
SO₄²⁻ (mM)
Ca²⁺ (mM)
Mg²⁺ (mM)
PO₄³⁻ (μM)
Ba²⁺ (μM)
DIC (mM)
δ¹³CDIC (‰, VPDB)

CL30

    55
0.0016
26.7
10.9
47.9
27.0
23.0
6.0
-20.4

    110
0.0012
24.4
10.5
48.3
31.5
19.8
7.6
-26.0

    170
0.0010
21.7
9.5
47.6
41.0
17.9
10.1
-31.9

    230
0.0051
19.0
8.9
47.5
49.9
17.6
11.1
-35.4

    290
0.0012
15.4
7.9
46.5
39.2
9.9
12.4
-38.0

    350
0.0000
13.1
7.0
46.1
45.5
20.0
14.3
-40.5

    410
0.0056
5.3
4.5
44.3
97.2
16.1
19.2
-45.7

    470
0.8752
1.0
3.7
43.8
112
46.5
22.8
-46.4

    530
3.2660
0.4
3.5
43.1
122
53.5
23.0
-40.0

    590
0.1318
0.2
3.4
43.7
116
60.8
23.7
-35.7

CL44

    10
0.0012
27.9
7.4
52.7
11.9
0.3
2.2
-10.9

    30
0.0004
27.9
7.4
52.4
10.1
0.3
2.6
-10.5

    50
0.0008
27.6
7.5
54.1
10.4
0.2
3.0
-12.1

    70
0.0003
27.2
7.0
52.7
9.3
0.2
3.0
-14.8

    90
0.0020
27.3
7.3
53.8
9.6
0.2
2.4
-10.9

    110
0.0005
26.7
7.2
54.0
9.1
0.2
2.4
-11.6

    130
0.0005
27.4
7.2
54.5
8.3
0.2
2.4
-11.9

    150
0.0008
26.9
7.3
54.3
7.1
0.2
2.5
-14.1

    170
0.0007
27.5
7.5
54.3
8.1
0.2
3.0
-13.1

    190
0.0006
27.5
7.3
54.7
10.9
0.2
2.5
-12.1

    210
0.0005
25.6
7.3
54.6
13.2
0.2
3.8
-17.2

    230
0.0017
25.9
7.0
54.2
13.7
0.2
3.5
-18.5

    250
0.0009
25.3
6.9
52.5
17.7
0.2
4.6
-20.4

    270
0.0018
24.1
6.3
51.6
22.3
0.3
5.5
-24.0

    290
0.0018
22.6
6.2
52.6
22.6
0.4
5.9
-23.1

    310
0.0014
21.5
6.0
53.0
35.1
0.3
6.8
-25.0

    330
0.0011
19.0
5.5
51.6
35.6
0.3
7.7
-28.8

    350
0.0017
19.7
5.1
51.0
47.3
0.3
8.1
-24.3

    370
0.0020
17.9
4.9
50.5
49.8
0.3
8.8
-30.1

    390
0.0015
16.4
4.6
50.5
50.8
0.4
9.5
-31.2

    410
0.0020
16.0
4.2
49.7
54.9
0.6
10.2
-31.9

    430
0.0012
15.6
3.9
47.8
52.1
0.5
11.0
-33.4

    450
0.0021
15.4
3.7
48.4
56.5
0.5
11.6
-35.1

    470
0.0022
13.2
3.6
50.3
55.4
0.6
12.0
-35.7

    490
0.0019
12.8
3.1
48.0
63.6
0.7
13.1
-35.8

    510
0.0026
12.0
2.7
48.4
67.1
0.8
13.1
-39.9

    530
0.0027
10.0
2.4
47.5
81.9
1.0
15.5
-39.1

    550
0.0024
8.9
2.2
47.4
88.0
1.3
16.3
-39.3

    570
0.0034
5.6
2.0
45.5
96.7
2.5
18.2
-40.5

    590
0.0039
4.7
1.8
44.6
101
4.9
19.5
-41.0

    610
0.0028
3.7
1.7
44.5
103
11.0
18.9
-39.7

    630
0.0022
2.7
1.6
44.2
57.0
19.3
20.5
-40.5

    650
0.0442
2.3
1.7
43.3
111
22.2
19.9
-39.6

    670
0.5139
2.1
1.8
44.7
109
31.7
21.1
-37.1

    690
1.1022
0.9
1.4
44.5
107
37.9
20.4
-37.4

    710
2.5450
0.8
1.3
43.9
102
38.6
21.5
-36.1

    730
0.0086
0.9
1.4
44.6
105
37.4
20
-36.8

    750
1.1475
1.1
1.2
43.0
86.8
36.0
19.2
-33.3

CL47

    55
0.0008
24.9
9.2
33.1
18.3
21.3
6.7
-21.1

    110
0.0006

21.1

    170
0.0004
26.0
10.9
39.7
16.4
13.6
6.8
-22.9

    230
0.0008

28.2

    290
0.0008
21.9
10.0
40.2
13.8
14.8
10.0
-28.5

    350
0.0009

44.7

    410
0.0009
16.5
8.8
39.8
43.1
13.9
13.7
-32.1

    470
0.0006
12.7
7.8
40.4
66.4
16.9
15.9
-34.0

    530
0.0006
9.3
6.5
40.0
65.9
13.7
18.2
-34.7

    590
0.0007
5.4
4.9
39.7
66.4
25.9
20.6
-38.2

    650
0.0773
1.1
4.5
39.8
76.3
52.0
24.2
-38.8

    710
0.7947
1.2
4.6
40.6
74.5
58.5
23.5
-37.4

    770
0.7692
1.4

70.2

22.1
-32.1

Vertical profiles of CL30, CL44, and CL47 for porewater barium concentrations and sediment barium contents together with barium/aluminium (Ba/Al) ratios are shown in Figure 5. Dissolved Ba²⁺ concentrations display maxima of 60.8, 38.6, and 58.5 μM below the SMTZ, respectively, and decreased upward towards the SMTZ (Figure 5). Bulk sediment Ba concentrations range from 306 to 957 mg kg⁻¹ (Table S6) with averages of 461 mg kg⁻¹ for CL30, 502 mg kg⁻¹ for CL44, and 502 mg kg⁻¹ for CL47. High Ba concentrations of bulk sediments at each core occur over narrow depth intervals (0.3–0.8 m) above the present SMTZ (Figure 5). Peak Ba concentrations within these zones reach 957, 741, and 790 mg kg⁻¹ and appear at approximately 4.3, 5.9, and 6.3 mbsf at cores CL30, CL44, and CL47, respectively. The refractory amount of solid phase barium at these cores amounts to 530, 550, and 590 mg kg⁻¹, respectively, which is considered to represent the “background” levels of solid phase barium [16]. Ba contents were normalized to Al in order to account for variations in lithology. Depth intervals with Ba content higher than these “background” levels are referred to as “Ba fronts.” At each core examined, the Ba fronts exist within 1.5 m above the depth of current SMTZ (Figure 5). The distance between the peak Ba concentration and the depth of sulfate depletion is approximately 0.4 m at CL30, 1.1 m at CL44, and 0.5 m at CL47.

(a)

(b)

(c)

Figure 5: Concentration depth profiles of dissolved barium (Ba²⁺), sulfate (SO₄²⁻), and solid-phase total barium (Batotal), barium/aluminium ratios (Ba/Al), and diagenetic barium (shown as red peaks) for cores CL30 (a), CL44 (b), and CL47 (c). Blue bands mark the SMTZ. Pink dash lines indicate the background barium contents based on the distribution of barium content. Barium contents above background represent diagenetic barite enrichments (red polygons).

POC contents at all the sites did not follow a general downward trend with average contents as 0.97% for CL30, 1.05% for CL44, and 1.05% for CL47 (Figures 2–4 and Table S6). The sediments in the study cores are mainly composed of silt and clay. At sites CL30 and CL44, the relative fractions of silt and clay are nearly constant with depth and the fractions of sand remain low values with depth except at ~320 cm in CL30 displaying elevated sand fraction (Figure S2). At site CL47, the sand fractions are low with depth in the interval of 0–200 cm, followed by an increase in sand fraction with two peaks at the depth of ~270 and ~430 cm. Below 500 cm, the sand fraction almost decreased to zero (Figure S2).

4.2. Timing of Authigenic Barite Front Accumulation

Dissolved barium fluxes towards the SMTZ were 1.58 mmol m⁻² yr⁻¹ for CL30, 1.54 mmol m⁻² yr⁻¹ for CL44, and 1.61 mmol m⁻² yr⁻¹ for CL47. The calculated time required for the formation of barite front is about 3.2, 3.0, and 1.3 kyr for the three cores, respectively, using an average porosity of 0.69 taken from CL44 (Table S5). Variations of porosity from 0.65 to 0.75 yield the time for barite front formation ranging between 2.2 and 4.2 kyr for CL30 and between 0.8 and 1.6 kyr for CL47. Sensitivity tests of the background Ba content, Ba²⁺ fluxes, and porosity are shown in Figures S5&S6.

4.3. Reaction-Transport Modeling

The modeled profiles and reaction rates are shown in Figures 2–4 and Table 4, respectively. The model parameters used to derive these results are listed in Tables S2-S4. The steady-state modeling reproduced the measured concentrations of SO₄²⁻, DIC, Ca²⁺, Mg²⁺, and PO₄³⁻ at sites CL44, CL47, and CL30 above the kink with obvious discrepancies between modeled and measured concentrations of CH₄ due to aforementioned degassing during core recovery (Figures 2–4). At site CL30, the model failed to reproduce the concentration gradients of SO₄²⁻, DIC, Ca²⁺, Mg²⁺, and PO₄³⁻ below the kink (~3.5 mbsf) which is likely caused by a transient condition that is not considered in the steady-state model.

Table 4: Depth-integrated simulated turnover rates and benthic methane fluxes based on the steady-state modeling.

CL30
CL44
CL47
Unit

: total POC mineralization rate
18.8
55.2
58.1
mmol m⁻² yr⁻¹ of C

: sulfate reduction via POC degradation
7.6
23.9
25.1
mmol m⁻² yr⁻¹ of SO₄²⁻

: methane formation via POC degradation
3.6
3.7
4.0
mmol m⁻² yr⁻¹ of CH₄

: gas dissolution
28.4
73.3
84.7
mmol m⁻² yr⁻¹ of CH₄

: anaerobic oxidation of methane
30.1
74.3
85.0
mmol m⁻² yr⁻¹ of CH₄

: authigenic CaCO₃ precipitation
3.0
2.2
5.6
mmol m⁻² yr⁻¹ of C

: authigenic MgCO₃ precipitation
5.4
7.1
0
mmol m⁻² yr⁻¹ of C

Sulfate consumed by AOM
79.8
75.7
77.2
%

Benthic flux of CH₄ at SWI
0.5
1.9
2.7
mmol m⁻² yr⁻¹ of CH₄

Percentage of CH₄ flux from depth
88.0
95.0
95.3
%

Percentage of CH₄ consumed by AOM
94.1
96.5
95.8
%

The sulfate concentration profile with a kink at site CL30 (Figure 2) could be explained by a recent increase in upward methane flux [9]. The linear extrapolation of the sulfate concentrations in the upper 3.5 m to zero sulfate concentration was taken as the initial condition for the non-steady-state model. Under this condition, the sulfate profile was fitted by a fixed CH₄ concentration (67 mM) at the lower boundary in equilibrium with the gas hydrate solubility under the conditions of in situ , , and . A sudden increase in CH₄ concentration reproduces the observed SO₄²⁻ concentration profile after running the model for ~85 yr (Figure 6). The increase in methane flux resulted in a prominent increase in the depth-integrated AOM rate from 30.1 mmol m⁻² y⁻¹ () to 140 mmol m⁻² y⁻¹ ().

Figure 6: Evolution of the sulfate profile over time from simulation of non-steady-state porewater profiles of core CL30.

The initial age of the organic matter was tuned until a good fit was obtained for the PO₄³⁻. The mean total depth-integrated rates of POC degradation were about 3 times higher at sites CL44 and CL47 (55.2 and 58.1 mmol m⁻² yr⁻¹) than that at site CL30 (18.8 mmol m⁻² yr⁻¹) (Table 4). The rates of POC degradation through sulfate reduction (POCSR) were 7.6, 23.9, and 25.1 mmol m⁻² yr⁻¹ at cores CL30, CL44, and CL47, respectively. In contrast to the relatively low rates of POCSR, AOM dominated the sulfate consumption with rates of 30.1, 74.3, and 84.7 mmol m⁻² yr⁻¹ for CL30, CL44, and CL47, respectively. The AOM rates were mainly sustained by an external methane source, and methanogenesis contributed only a negligible amount of methane (Table 4). The AOM consumes almost all the CH₄ with benthic CH₄ fluxes of 0.49, 2.0, and 2.7 mmol m⁻² yr⁻¹ at sites CL30, CL44, and CL47, respectively.

5. Discussion

5.1. Formation Mechanisms of the Nonlinear Porewater Profiles

The sulfate concentration profiles of the porewater in marine sediments depend on the availability of labile organic matter amenable to sulfate reducers, diffusive/advective methane flux, and depositional conditions [9, 59–64]. Combination of these factors would result in linear, kink, concave-up, concave-down, and sigmoidal (S-shape) type sulfate concentration trends in marine sediments [9].

The porewater profiles of the three study cores exhibit kink-type features. The plausible mechanisms for the occurrence of the kink-type profile include (1) irrigation and seawater intrusion due to biological, physical, and hydrological processes; (2) changes in the sedimentation rate or porosity due to depositional events; and (3) changes in methane flux and upward advection of fluid [14]. Bioirrigation has been shown to generally occur in a decimeter scale in the surface sediments [65, 66]. In fact, no macroorganisms were observed in the study cores below the upper few centimeters of sediment. The lithology of the upper two-three meters of the sediments in the study cores was dominated by fine-grained hemipelagic sediments mainly consisting of silty clay without any discernible abnormal deposition (Figure S2). Although deep-water turbidity current channel and fan systems are well developed in the study region [67], the homogeneous grain size distributions in cores CL44 and CL47 reveal that the sediment above the kinks of sulfate was not impacted by turbidites, which are typically characterized by upward grading in grain size. The C-M plot also suggests the absence of turbidites in the study cores (Figure S2 [68]). Moreover, by comparing the depth profiles of CaCO₃ content in CL44 and CL47 with that in an adjacent core (SO49-37KL) with established Marine Isotope Stage, we found that the upper ~2 m sediments of CL44 and CL47 represent normal hemipelagic background deposition during the Holocene (Figure S3) [69, 70]. The relatively constant ratios of Ti/Al, Si/Al, and Zr/Rb above the kinks indicate a stable input of detrital fraction (Figure S4). In contrast, the layers at the interval of ~1.4 to ~4.2 mbsf in CL44 and ~1.8 to ~5 mbsf in CL47 exhibiting high Si, Ti, and Zr/Rb contents and coarser grain sizes (Figure S4) suggest elevated input of detrital fraction during sea-level lowstands [71, 72]. In addition, the flat seafloor topography in the study area also precludes the occurrence of abrupt depositional event such as landslide (Figure 1). Therefore, it is unlikely that the irrigation-like feature in CL44 and CL47 was caused by mass-transport deposits [44]. Furthermore, there is no indicator for upward fluid advection at sites CL44 and CL47.

We argue that the cause for the formation of the irrigation-like porewater profiles is probably the bubble irrigation by rising free gas through escaping tubes [7, 12, 51]. Such features were observed at the nearby Haima cold seeps and attributed to bubble irrigation or a recent increase in methane flux [33]. Moreover, BSR and acoustic blanking which are indicative of free gas accumulation were identified in the study area (Fang Y., unpublished data). Hence, gas bubble irrigation is the most likely mechanism to explain the observed profiles at cores CL44 and CL47.

At core CL30, the sediments consist of homogenous silty clay without discernible abnormal deposition and the sulfate concentrations decrease gradually without maintaining seawater-like values above the kink at 3.5 mbsf (Figure 2). We thus hypothesize that the kink in the sulfate profile at core CL30 results from a (sub)recent increase in the upward methane flux, similar to the scenario reported in the Sea of Marmara, the continental margin offshore Pakistan, the slope area south of Svalbard, the Niger Delta, the southern SCS, and so on [10–13, 44]. A simplified numerical model exercise, assuming a diffusional porewater system with POCSR and AOM as the only biogeochemical reactions, was used to demonstrate this scenario (Figure 6). The assumption of diffusive transport of porewater species is warranted because it has been suggested that porewater solute distributions are dominated by diffusion even if free gas transport and fluid advection exist [14, 73].

The current barite front is located at about 4.2-6.4 mbsf, very close to the current SMTZ (4.7-7 mbsf), indicating that the barite front might form in the recent past to the present day induced by a recent enhancement of methane flux [11]. Actually, the measured SO₄²⁻ concentration profile can be reproduced after a sudden increase in CH₄ concentration lasting for ~85 yr. On the other hand, based on the calculated diffusive Ba²⁺ fluxes and the depth-integrated Ba contents, the time required to form the observed authigenic barite front above the current SMTZ is about 2.2-4.2 kyr for CL30, given the uncertainties of porosity. The difference of estimated duration of constant methane flux between these two approaches may suggest that the barite front was not a result of the recent increase in methane flux inducing the kink-type sulfate profile. Instead, it is more likely that the SMTZ has experienced several fluctuations in depth, considering the episodic pulses of upward methane flux which have occurred in this area as shown by previous studies [25, 32, 33]. However, this decoupled record between sediments and porewaters is commonly observed at cold seeps [74–77] and is considered to reflect the variations of methane fluxes and the resulting SMTZ in the sedimentary column [74]. Observations and numerical modeling suggest that the response of porewater geochemical signatures is more rapid on timescales of months to centuries than the accumulation of authigenic barite deposits on timescales of decades to hundred thousands of years [11, 12, 14, 16–19, 74]. On the whole, our results suggest that combining porewater data with sedimentary barite front records may provide important clues for better understanding of the evolution of methane seepage.

5.2. Methane-Related Carbon Cycling and Source of Methane

Based on the simulation results derived from the steady-state modeling, AOM consumed ~80%, 76%, and 77% of sulfate in CL30, CL44, and CL47, respectively. AOM thus acts as an efficient barrier preventing methane from being released into the water column at the studied cores. This is supported by the low δ¹³CDIC values at the SMTZs mainly derived from methane. AOM increases porewater alkalinity by producing bicarbonate and results in the precipitation of authigenic carbonates as shown by the decrease in Ca²⁺ and Mg²⁺ concentrations with depth (Figures 3–6).

In addition, the δ¹³CDIC values below the SMTZ become more positive than those at the SMTZ. The reversal in δ¹³CDIC below the SMTZ is caused by the generation of ¹³C-enriched DIC via local methanogenesis at the methanogenic zone [63, 78]. The ¹³C-enriched DIC would migrate into the SMTZ from the methanogenic zone and “dilute” the ¹²C pool of DIC in porewater. Thus, in a closed system, DIC generated by local methanogenesis is an important source of DIC in the carbon budget within the SMTZ [78, 79].

Based on the modeling results of methane turnovers, the depth-intergrated AOM rates at cores CL30, CL44, and CL47 are about 8 to 21 times of the in situ methanogenesis rates (Table 4). Therefore, the relative proportions of external methane sources contributed to the total methane pool are 88%, 95%, and 95% at cores CL30, CL44, and CL47, respectively. This indicates that the majority of methane fuelling AOM at the SMTZ was sourced from subsurface sediments. There are two general pathways for producing methane in marine sediments, including microbial methane generated via CO₂ reduction or the fermentation of reduced carbon substrates (e.g., acetate and methanol; [80]) and thermogenic methane formed via thermal cracking of organic matter and/or heavy hydrocarbons [81]. The δ¹³C values of methane are generally distinct between these two types of methane. The δ¹³C values of microbial methane typically range from −50‰ to −110‰ [80], whereas those of thermogenic methane range from −30‰ to −50‰ [81].

Because δ¹³C values of headspace methane in the sediments are absent in the study area, porewater DIC content and δ¹³CDIC are utilized to constrain the origin of methane. Generally, porewater DIC in marine sediments is mainly derived from (1) the DIC that is diffusing from the overlying seawater into the sediments or the seawater DIC trapped within sediments during burial, (2) the DIC generated by the degradation of sedimentary organic matter, (3) the DIC produced by AOM, and (4) the residual DIC derived from methanogenesis [82, 83]. In order to obtain the carbon isotopic composition of DIC derived from external methane, we applied a simple four-end-member mixing model. The four end-members are (1) seawater-derived DIC trapped within sediments during burial (SW), (2) DIC produced by POCSR, (3) DIC derived from external methane (EM) via AOM, and (4) DIC generated by in situ methanogenesis (ME). Note that methane production via local methanogenesis was assumed to be competently recycled by AOM. As a result, the carbon isotopic composition of DIC produced by local methanogenesis was identical to that of organic matter (OM) [82–84]. In a closed system, δ¹³C balance of the porewater DIC pool at the SMTZ can be expressed by

where is the proportion of DIC contributed to the total DIC pool and the subscripts , , , and refer to DIC derived from seawater, organic matter, external methane, and local methanogenesis, respectively. The values are estimated as typical seawater DIC concentration (2.1 mM) divided by DIC concentration at the SMTZ, and the δ¹³CSW is assumed to be 0‰. The δ¹³C value of sedimentary organic matter in sediment of the SCS (−20‰; [85]) is used for the δ¹³COM. The overall δ¹³C of DIC derived from methanogenesis is equal to the δ¹³COM assuming that methane produced by local methanogenesis was completely converted to DIC by AOM [84]. The contribution fractions of OSR, AOM, and local methanogenesis shown as , , and are calculated from the steady-state modeling. The δ¹³Cex can be acquired from a regression of porewater vs. DIC (Figure 7). The regression commonly shows linear in seep-impacted sediments, thus providing definitive δ¹³CDIC supplied to porewater [82, 86, 87]. The δ¹³Cex values and the contribution fractions estimated from the model are listed in Table 5. The estimated δ¹³CEM values in the shallow sediments are -74.1‰ (CL30), -75.4‰ (CL44), and -66.7‰ (CL47), suggesting the external methane migrating into the shallow sediments is microbial in origin [83, 88]. The absence of higher hydrocarbons in headspace gas samples also supports the microbial origin of methane in the study area.

(a)

(b)

(c)

Figure 7: Plots of DIC vs. . The δ¹³Cex were calculated using the linear regression of vs. DIC for CL30 (a), CL44 (b), and CL47 (c).

Table 5: Fractions of DIC from different sources contributed to the total DIC pool and δ¹³C values from different sources.

Core ID
δ¹³CSMTZ (‰, VPDB)
δ¹³Cex (‰, VPDB)
δ¹³CSW (‰, VPDB)
δ¹³COM (‰, VPDB)
(%)
(%)
(%)
(%)
δ¹³Cmethane (‰, VPDB)

CL30
-46.4
-48.1
0
-20
9.8
6.6
28.0
55.5
-74.1

CL44
-41.0
-46.1
0
-20
9.8
5.3
34.3
50.6
-75.4

CL47
-38.8
-43.1
0
-20
8.7
5.2
32.9
53.1
-66.7

Previous studies have suggested that microbial methane was the main hydrocarbon source with minor contributions of oil-derived compounds and pyrolysis gas at the Haima cold seeps [32, 33, 89, 90]. Gas chimney structures, which are well developed around the Haima cold seeps, might serve as conduits for the upward migration of biogenic gas from the underlying free gas reservoir beneath the GHSZ to shallow sediments [29, 42]. This observation suggests that microbial methane at the study sites might be derived from an underlying free gas reservoir trapped beneath the GHSZ.

5.3. Implications of the Time Constraint on Methane Seepage

Based on the diffusive dissolved barium flux and the excess barium content in the sediments, the time for the observed barite enrichments just above the current SMTZ is estimated to be about 3 kyr and 0.8-1.6 kyr for CL44 and CL47, respectively, given the uncertainties of porosity. These results suggest that the SMTZ has been fixed at the current sediment depth for a time period of at least several thousand years at these sites. The irrigation-type sulfate profiles are possibly maintained by continuous mixing of seawater into the sediment over these time periods, like the case in the pockmark sediments of Congo Fan [74]. Furthermore, the depth of SMTZ was speculated to have fluctuated due to variations in methane flux as suggested by the difference in the estimated duration of the barite enrichment and the recent increase methane flux at CL30. Overall, our results show that the methane flux has been fluctuating over the last hundreds to thousands of years in the vicinity of Haima cold seeps.

In fact, methane seepages around the Haima cold seeps are characterized by distinct periodicity of seep activities during the past several thousands of years. Radiocarbon ages of bivalve shells suggest that a major seepage event occurred during the period of 6.1 to 5.1 ka B.P., followed by a subordinate seepage event spanning from 3.9 to 2.9 ka B.P. at the Haima cold seeps [25]. The widespread occurrence of dead bivalves on the seafloor reflects a decline in current seepage intensity [25]. Moreover, modeling of porewater profiles at the Haima cold seeps predicts that gas hydrate formation in the seepage center started at least 150 yr B.P. and the subsequent sealing of gas hydrates favored the lateral migration of methane-rich fluids in the coarser, more permeable interval [33]. Sedimentation dynamics, including sediment instabilities and mass wasting, may trigger the destabilization of the gas hydrate reservoir and the resulting occurrence of methane seepage. The evolution and fate of methane seepage are also considered to be affected by local fluid flow dynamics and associated migration of both free gas and methane-rich fluids along fractures, as well as the redirection of gas supply from the reservoir due to pore space clogging by gas hydrate in shallow sediments [25, 33]. The exact mechanism of the changes in methane flux around the Haima cold seeps area is beyond the scope of this study. Despite this, our quantitative study provides some constraint on the duration of methane seepage and may have implication for understanding the evolution of methane seepage in the petroliferous Qiongdongnan Basin.

6. Conclusions

This study is aimed at understanding the methane source and turnover as well as provide some constraints on the timing of methane seepage to the west of “Haima cold seeps.” The steady-state reaction-transport modeling of SO₄²⁻, CH₄, DIC, PO₄³⁻, Ca²⁺, and Mg²⁺ in CL44 and CL47 suggests that gas bubble transport may lead to the irrigation-like feature in the upper 2 m and relatively high AOM rates (74.3 mmol m⁻² yr⁻¹ for CL44 and 85.0 mmol m⁻² yr⁻¹ for CL47). The time required for the enrichment of authigenic barium fronts slightly above the current SMTZ is approximately 3 kyr for CL44 and 0.8-1.6 kyr for CL47, respectively. In contrast, a recent increase in methane flux (prior to ~85 yr) is the likely cause of the kink at 3.5 m of the sulfate profile in CL30 demonstrated by the transient-state modeling. The estimated time required for the formation of the diagenetic barium peak just above the current SMTZ was 2.2-4.2 kyr at this core. The discrepancy in the time estimates constrained by two different approaches suggests that the position of SMTZ possibly has fluctuated due to variation in methane flux at the site. In addition, based on the four DIC end-member mixing calculation, the δ¹³C values of the external methane in cores CL30, CL44, and CL47 are -74.1‰, -75.4‰, and -66.7‰, respectively. This is indicative of the biogenic origin of external methane from an underlying reservoir. Our results suggest that methane seepage exists in a broader area in the vicinity of the “Haima cold seeps” and the methane fluxes may have fluctuated frequently for the last several hundreds to thousands of years.

Data Availability

The data used to support the findings of this study are included within the main text and the supplementary materials.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

We thank the crew of the Haiyang-4 exploration ship for collecting the piston cores. We are also grateful to our colleagues from GMGS for their collection and analysis to the porewater and sediment samples. The authors appreciate Dr. T. Yang (Nanjing University) and Dr. X. J. Yin (Third Institute of Oceanography, Ministry of Natural Resources) for the help with the geochemical analyses. This study was supported by the National Key R&D Program of China (Grant: 2018YFC0310006), the China Postdoctoral Science Foundation (Grant: 2017M622654), the open-funds of Key Laboratory of Marine Mineral Resources, Ministry of Land and Resources (KLMMR-2017-A-08), the National Special Project on Gas Hydrate of China (Grant: GZH201100301), and the National Natural Science Foundation of China (Grant: 41730528, 41806074).

Supplementary Materials

In Table S1, reaction terms of all species used in the model are listed. Table S2 presents the summary of model parameters and boundary conditions used in the steady-state model, whereas Table S3 lists the derived and measured parameters used in the steady-state model. Parameters used in the nonsteady-state model are listed in Table S4. In Table S5, porosity and dry grain density of core CL44 are listed. Table S6 lists the solid-phase element contents, Ba/Al, Ti/Al, Si/Al, and Zr/Rb ratios, as well as POC and CaCO₃ contents of cores CL30, CL44, and CL47. Table S7 shows the grain size parameters of the studied cores. In addition, Figure S1 displays the measured (symbols) and modelled (curve) porosity of core CL44. Figure S2 shows the grain size parameters of the sediment cores. Figure S3 displays the downcore variations in carbonate contents and the age models of the studied cores. Figure S4 shows the Ti/Al, Si/Al, Zr/Rb, and carbonate contents of the sediments cores. Figure S5 details the uncertainty evaluation of the formation times of the diagenetic barium enrichments considering variable solid-phase barium background contents and porosities, whereas Figure S6 shows the uncertainty evaluation of the formation times of the diagenetic barium enrichments considering variable Ba²⁺ fluxes and porosities. In summary, Tables S1–S4 are the detailed descriptions of the reaction-transport models. Tables S5–S7 and Figure S1 are the datasets of the geotechnical and geochemical proxies of sediment samples, which contribute to the concentration depth profiles in Figure 5 of this manuscript. Figures S2–S4 show downcore variations in the geotechnical and geochemical proxies of sediment cores. Figures S5 and S6 are the sensitivity analysis of the formation times of the diagenetic barium enrichments. (Supplementary Materials)

References

A. G. Judd and M. Hovland, Submarine Fluid Flow, the Impact on Geology, Biology, and the Marine Environment, Cambridge University Press, Cambridge, 2007.

T. Yuan, G. D. Spence, R. D. Hyndman, T. A. Minshull, and S. C. Singh, “Seismic velocity studies of a gas hydrate bottom-simulating reflector on the northern Cascadia continental margin: amplitude modeling and full waveform inversion,” Journal of Geophysical Research: Solid Earth, vol. 104, no. B1, pp. 1179–1191, 1999.

R. O. Barnes and E. D. Goldberg, “Methane production and consumption in anoxic marine sediments,” Geology, vol. 4, no. 5, pp. 297–300, 1976.

A. Boetius, K. Ravenschlag, C. J. Schubert et al., “A marine microbial consortium apparently mediating anaerobic oxidation of methane,” Nature, vol. 407, no. 6804, pp. 623–626, 2000.

R. Luff, J. Greinert, K. Wallmann, I. Klaucke, and E. Suess, “Simulation of long-term feedbacks from authigenic carbonate crust formation at cold vent sites,” Chemical Geology, vol. 216, no. 1-2, pp. 157–174, 2005.

W. Ussler III and C. K. Paull, “Rates of anaerobic oxidation of methane and authigenic carbonate mineralization in methane-rich deep-sea sediments inferred from models and geochemical profiles,” Earth and Planetary Science Letters, vol. 266, no. 3-4, pp. 271–287, 2008.

M. Haeckel, B. P. Boudreau, and K. Wallmann, “Bubble-induced porewater mixing: a 3-D model for deep porewater irrigation,” Geochimica et Cosmochimica Acta, vol. 71, no. 21, pp. 5135–5154, 2007.

S. Sommer, O. Pfannkuche, P. Linke et al., “Efficiency of the benthic filter: biological control of the emission of dissolved methane from sediments containing shallow gas hydrates at Hydrate Ridge,” Global Biogeochemical Cycles, vol. 20, no. 2, 2006.

C. Hensen, M. Zabel, K. Pfeifer et al., “Control of sulfate pore-water profiles by sedimentary events and the significance of anaerobic oxidation of methane for the burial of sulfur in marine sediments,” Geochimica et Cosmochimica Acta, vol. 67, no. 14, pp. 2631–2647, 2003.

P. Halbach, E. Holzbecher, T. Reichel, and R. Moche, “Migration of the sulphate–methane reaction zone in marine sediments of the Sea of Marmara—can this mechanism be tectonically induced?” Chemical Geology, vol. 205, no. 1-2, pp. 73–82, 2004.

D. Fischer, J. M. Mogollón, M. Strasser et al., “Subduction zone earthquake as potential trigger of submarine hydrocarbon seepage,” Nature Geoscience, vol. 6, no. 8, pp. 647–651, 2013.

W. L. Hong, M. E. Torres, J. L. Carroll et al., “Seepage from an arctic shallow marine gas hydrate reservoir is insensitive to momentary ocean warming,” Nature Communications, vol. 8, no. 1, article 15745, 2017.

A. de Prunelé, L. Ruffine, V. Riboulot et al., “Focused hydrocarbon-migration in shallow sediments of a pockmark cluster in the Niger Delta (Off Nigeria),” Geochemistry, Geophysics, Geosystems, vol. 18, no. 1, pp. 93–112, 2017.

W. L. Hong, S. Sauer, G. Panieri et al., “Removal of methane through hydrological, microbial, and geochemical processes in the shallow sediments of pockmarks along eastern Vestnesa Ridge (Svalbard),” Limnology and Oceanography, vol. 61, no. S1, pp. S324–S343, 2016.

M. E. Torres, G. Bohrmann, and E. Suess, “Authigenic barites and fluxes of barium associated with fluid seeps in the Peru subduction zone,” Earth and Planetary Science Letters, vol. 144, no. 3-4, pp. 469–481, 1996.

G. R. Dickens, “Sulfate profiles and barium fronts in sediment on the Blake Ridge: present and past methane fluxes through a large gas hydrate reservoir,” Geochimica et Cosmochimica Acta, vol. 65, no. 4, pp. 529–543, 2001.

N. Riedinger, S. Kasten, J. Gröger, C. Franke, and K. Pfeifer, “Active and buried authigenic barite fronts in sediments from the Eastern Cape Basin,” Earth and Planetary Science Letters, vol. 241, no. 3-4, pp. 876–887, 2006.

G. T. Snyder, A. Hiruta, R. Matsumoto et al., “Pore water profiles and authigenic mineralization in shallow marine sediments above the methane-charged system on Umitaka Spur, Japan Sea,” Deep Sea Research Part II Topical Studies in Oceanography, vol. 54, no. 11-13, pp. 1216–1239, 2007.

S. Kasten, K. Nöthen, C. Hensen, V. Spieß, M. Blumenberg, and R. R. Schneider, “Gas hydrate decomposition recorded by authigenic barite at pockmark sites of the northern Congo Fan,” Geo-Marine Letters, vol. 32, no. 5-6, pp. 515–524, 2012.

L. F. Fan, S. Lin, C. W. Hsu, Y. T. Tseng, T. F. Yang, and K. M. Huang, “Formation and preservation of authigenic pyrite in the methane dominated environment,” Deep-Sea Research Part I, vol. 138, pp. 60–71, 2018.

D. F. Chen, Y. Y. Huang, X. L. Yuan, and L. M. Cathles III, “Seep carbonates and preserved methane oxidizing archaea and sulfate reducing bacteria fossils suggest recent gas venting on the seafloor in the northeastern South China Sea,” Marine and Petroleum Geology, vol. 22, no. 5, pp. 613–621, 2005.

X. Han, E. Suess, Y. Huang et al., “Jiulong methane reef: microbial mediation of seep carbonates in the South China Sea,” Marine Geology, vol. 249, no. 3-4, pp. 243–256, 2008.

H. Tong, D. Feng, H. Cheng et al., “Authigenic carbonates from seeps on the northern continental slope of the South China Sea: new insights into fluid sources and geochronology,” Marine and Petroleum Geology, vol. 43, pp. 260–271, 2013.

D. Feng and D. Chen, “Authigenic carbonates from an active cold seep of the northern South China Sea: new insights into fluid sources and past seepage activity,” Deep Sea Research Part II Topical Studies in Oceanography, vol. 122, pp. 74–83, 2015.

Q. Liang, Y. Hu, D. Feng et al., “Authigenic carbonates from newly discovered active cold seeps on the northwestern slope of the South China Sea: constraints on fluid sources, formation environments, and seepage dynamics,” Deep Sea Research Part I Oceanographic Research Papers, vol. 124, pp. 31–41, 2017.

D. Feng, J. W. Qiu, Y. Hu et al., “Cold seep systems in the South China Sea: an overview,” Journal of Asian Earth Sciences, vol. 168, pp. 3–16, 2018.

Y. Bai, H. Song, Y. Guan, J. Chen, and B. Liu, “Structural characteristics and genesis of pockmarks in the northwest of the South China Sea derived from reflective seismic and multibeam data,” Chinese Journal of Geophysics, vol. 57, pp. 2208–2222, 2014.

B. Liu and S. Liu, “Gas bubble plumes observed at north slope of South China Sea from multi-beam water column data,” Haiyang Xuebao, vol. 39, pp. 83–89, 2017.

L. Yang, B. Liu, M. Xu, S. Liu, Y. Guan, and Y. Gu, “Characteristics of active cold seepages in Qiongdongnan Sea Area of the northern South China Sea,” Chinese Journal of Geophysics, vol. 61, pp. 2905–2914, 2018.

B. Zhao, S. Liu, L. Li, and J. Guo, “Distribution pattern of cold seeps in South China Sea and its geological significance,” Marine Geology Frontiers, vol. 34, pp. 32–43, 2018.

J. Ye, J. Wei, J. Liang et al., “Complex gas hydrate system in a gas chimney, South China Sea,” Marine and Petroleum Geology, vol. 104, pp. 29–39, 2019.

X. Wang, N. Li, D. Feng et al., “Using chemical compositions of sediments to constrain methane seepage dynamics: a case study from Haima cold seeps of the South China Sea,” Journal of Asian Earth Sciences, vol. 168, pp. 137–144, 2018.

Y. Hu, M. Luo, Q. Liang et al., “Pore fluid compositions and inferred fluid flow patterns at the Haima cold seeps of the South China Sea,” Marine and Petroleum Geology, vol. 103, pp. 29–40, 2019.

B. Taylor and D. E. Hayes, “The tectonic evolution of the South China Basin,” in The Tectonic and Geologic Evolution of Southeast Asian Seas and Islands, D. E. Hayes, Ed., vol. 23 of Geophysical Monograph Series, pp. 89–104, AGU, 1980.

K. Ru and J. D. Pigott, “Episodic rifting and subsidence in the South China Sea,” AAPG Bulletin, vol. 70, pp. 1136–1155, 1986.

T. Lüdmann and H. K. Wong, “Neotectonic regime on the passive continental margin of the northern South China Sea,” Tectonophysics, vol. 311, no. 1-4, pp. 113–138, 1999.

X. Xie, R. D. Müller, S. Li, Z. Gong, and B. Steinberger, “Origin of anomalous subsidence along the northern South China Sea margin and its relationship to dynamic topography,” Marine and Petroleum Geology, vol. 23, no. 7, pp. 745–765, 2006.

X. Shi, H. Jiang, J. Yang, X. Yang, and H. Xu, “Models of the rapid post-rift subsidence in the eastern Qiongdongnan Basin, South China Sea: implications for the development of the deep thermal anomaly,” Basin Research, vol. 29, no. 3, pp. 340–362, 2017.

Q. Sun, J. Cartwright, T. Lüdmann, S. Wu, and G. Yao, “Three-dimensional seismic characterization of a complex sediment drift in the South China Sea: evidence for unsteady flow regime,” Sedimentology, vol. 64, no. 3, pp. 832–853, 2017.

X. Wang, S. Wu, S. Yuan et al., “Geophysical signatures associated with fluid flow and gas hydrate occurrence in a tectonically quiescent sequence, Qiongdongnan Basin, South China Sea,” Geofluids, vol. 10, no. 3, 368 pages, 2010.

G. Hui, S. Li, L. Guo et al., “Source and accumulation of gas hydrate in the northern margin of the South China Sea,” Marine and Petroleum Geology, vol. 69, pp. 127–145, 2016.

J. Wang, S. Wu, X. Kong et al., “Subsurface fluid flow at an active cold seep area in the Qiongdongnan Basin, northern South China Sea,” Journal of Asian Earth Sciences, vol. 168, pp. 17–26, 2018.

K. Grasshoff, K. Kremling, and M. Ehrhardt, “Methods of Seawater Analysis,” Third, Completely Revised and Extended Edition, Wiley-VCH, Weinheim, Germany, 1999.

J. Feng, S. Yang, J. Liang et al., “Methane seepage inferred from the porewater geochemistry of shallow sediments in the Beikang Basin of the southern South China Sea,” Journal of Asian Earth Sciences, vol. 168, pp. 77–86, 2018.

T. Yang, S. Y. Jiang, J. H. Yang et al., “Dissolved inorganic carbon (DIC) and its carbon isotopic composition in sediment pore waters from the Shenhu area, northern South China Sea,” Journal of Oceanography, vol. 64, no. 2, pp. 303–310, 2008.

H. D. Schulz, “Quantification of early diagenesis: dissolved constituents in pore water and signals in the solid phase,” in Marine Geochemistry, H. D. Schulz and M. Zabel, Eds., pp. 73–124, Springer, Berlin, Germany, 2006.

B. P. Boudreau, Diagenetic Models and Their Implementation: Modelling Transport and Reactions in Aquatic Sediments, Springer-Verlag, Berlin, 1997.

K. Wallmann, G. Aloisi, M. Haeckel, A. Obzhirov, G. Pavlova, and P. Tishchenko, “Kinetics of organic matter degradation, microbial methane generation, and gas hydrate formation in anoxic marine sediments,” Geochimica et Cosmochimica Acta, vol. 70, no. 15, pp. 3905–3927, 2006.

P. C. Chuang, A. W. Dale, K. Wallmann et al., “Relating sulfate and methane dynamics to geology: accretionary prism offshore SW Taiwan,” Geochemistry, Geophysics, Geosystems, vol. 14, no. 7, pp. 2523–2545, 2013.

M. Luo, A. W. Dale, K. Wallmann et al., “Estimating the time of pockmark formation in the SW Xisha Uplift (South China Sea) using reaction-transport modeling,” Marine Geology, vol. 364, pp. 21–31, 2015.

M. Luo, A. W. Dale, L. Haffert et al., “A quantitative assessment of methane cycling in Hikurangi Margin sediments (New Zealand) using geophysical imaging and biogeochemical modeling,” Geochemistry, Geophysics, Geosystems, vol. 17, no. 12, pp. 4817–4835, 2016.

Z. Duan, N. Møller, J. Greenberg, and J. H. Weare, “The prediction of methane solubility in natural waters to high ionic strength from 0 to 250 °C and from 0 to 1600 bar,” Geochimica et Cosmochimica Acta, vol. 56, no. 4, pp. 1451–1460, 1992.

J. J. Middelburg, “A simple rate model for organic matter decomposition in marine sediments,” Geochimica et Cosmochimica Acta, vol. 53, no. 7, pp. 1577–1581, 1989.

M. J. Whiticar, E. Faber, and M. Schoell, “Biogenic methane formation in marine and freshwater environments: CO₂ reduction vs. acetate fermentation–isotope evidence,” Geochimica et Cosmochimica Acta, vol. 50, no. 5, pp. 693–709, 1986.

E. A. Solomon, A. J. Spivack, M. Kastner, M. E. Torres, and G. Robertson, “Gas hydrate distribution and carbon sequestration through coupled microbial methanogenesis and silicate weathering in the Krishna-Godavari Basin, offshore India,” Marine and Petroleum Geology, vol. 58, pp. 233–253, 2014.

F. J. Millero, “Thermodynamics of the carbon dioxide system in the oceans,” Geochimica et Cosmochimica Acta, vol. 59, no. 4, pp. 661–677, 1995.

R. E. Zeebe and D. A. Wolf-Gladrow, CO2 in Seawater: Equilibrium, Kinetics and Isotopes, Elsevier, U. K., 2001.

P. Tishchenko, C. Hensen, K. Wallmann, and C. S. Wong, “Calculation of the stability and solubility of methane hydrate in seawater,” Chemical Geology, vol. 219, no. 1-4, pp. 37–52, 2005.

C. Niewohner, C. Hensen, S. Kasten, M. Zabel, and H. D. Schulz, “Deep sulfate reduction completely mediated by anaerobic methane oxidation in sediments of the upwelling area off Namibia,” Geochimica et Cosmochimica Acta, vol. 62, no. 3, pp. 455–464, 1998.

M. Zabel and H. D. Schulz, “Importance of submarine landslides for non-steady state conditions in pore water systems — lower Zaire (Congo) deep-sea fan,” Marine Geology, vol. 176, no. 1-4, pp. 87–99, 2001.

M. Schmidt, C. Hensen, T. Mörz et al., “Methane hydrate accumulation in “Mound 11” mud volcano, Costa Rica forearc,” Marine Geology, vol. 216, no. 1-2, pp. 83–100, 2005.

T. Tréude, J. Niggemann, J. Kallmeyer et al., “Anaerobic oxidation of methane and sulfate reduction along the Chilean continental margin,” Geochimica et Cosmochimica Acta, vol. 69, no. 11, pp. 2767–2779, 2005.

M. Kastner, G. Claypool, and G. Robertson, “Geochemical constraints on the origin of the pore fluids and gas hydrate distribution at Atwater Valley and Keathley canyon, northern Gulf of Mexico,” Marine and Petroleum Geology, vol. 25, no. 9, pp. 860–872, 2008.

W. L. Hong, E. A. Solomon, and M. E. Torres, “A kinetic-model approach to quantify the effect of mass transport deposits on pore water profiles in the Krishna–Godavari Basin, Bay of Bengal,” Marine and Petroleum Geology, vol. 58, pp. 223–232, 2014.

R. C. Aller and J. Y. Aller, “Meiofauna and solute transport in marine muds,” Limnology and Oceanography, vol. 37, no. 5, pp. 1018–1033, 1992.

C. Meile, C. M. Koretsky, and P. V. Cappellen, “Quantifying bioirrigation in aquatic sediments: an inverse modeling approach,” Limnology and Oceanography, vol. 46, no. 1, pp. 164–177, 2001.

S. Yuan, S. Wu, T. Lüdmann et al., “Fine-grained Pleistocene deepwater turbidite channel system on the slope of Qiongdongnan Basin, northern South China Sea,” Marine and Petroleum Geology, vol. 26, no. 8, pp. 1441–1451, 2009.

R. Passega, “Grain size representation by CM patterns as a geologic tool,” Journal of Sedimentary Petrology, vol. 34, no. 4, pp. 830–847, 1964.

J. Schönfeld and H. R. Kudrass, “Hemipelagic sediment accumulation rates in the South China Sea related to late Quaternary sea-level changes,” Quaternary Research, vol. 40, no. 3, pp. 368–379, 1993.

J. Qian, Paleooceanoraphy for the Late Quaternary in the South China Sea, China Science Press, Beijing, 1999, (in Chinese).

Y. Hu, L. Chen, D. Feng, Q. Liang, Z. Xia, and D. Chen, “Geochemical record of methane seepage in authigenic carbonates and surrounding host sediments: a case study from the South China Sea,” Journal of Asian Earth Sciences, vol. 138, pp. 51–61, 2017.

Y. Xu, L. Liu, H. Zhou et al., “Submarine landslide identified in DLW3102 core of the northern continental slope, South China Sea,” Journal of Ocean University of China, vol. 17, no. 1, pp. 147–155, 2018.

J. H. Kim, M. E. Torres, J. Choi, J. J. Bahk, M. H. Park, and W. L. Hong, “Inferences on gas transport based on molecular and isotopic signatures of gases at acoustic chimneys and background sites in the Ulleung Basin,” Organic Geochemistry, vol. 43, pp. 26–38, 2012.

K. Nöthen and S. Kasten, “Reconstructing changes in seep activity by means of pore water and solid phase Sr/Ca and Mg/Ca ratios in pockmark sediments of the northern Congo Fan,” Marine Geology, vol. 287, no. 1-4, pp. 1–13, 2011.

A. Mazumdar, P. Dewangan, H. M. Joäo et al., “Evidence of paleo-cold seep activity from the Bay of Bengal, offshore India,” Geochemistry, Geophysics, Geosystems, vol. 10, no. 6, 2009.

Y. C. Lim, S. Lin, T. F. Yang, Y. G. Chen, and C. S. Liu, “Variations of methane induced pyrite formation in the accretionary wedge sediments offshore southwestern Taiwan,” Marine and Petroleum Geology, vol. 28, no. 10, pp. 1829–1837, 2011.

L. M. Wehrmann, S. P. Templer, B. Brunner, S. M. Bernasconi, L. Maignien, and T. G. Ferdelman, “The imprint of methane seepage on the geochemical record and early diagenetic processes in cold-water coral mounds on Pen Duick Escarpment, Gulf of Cadiz,” Marine Geology, vol. 282, no. 1-2, pp. 118–137, 2011.

S. Chatterjee, G. R. Dickens, G. Bhatnagar et al., “Pore water sulfate, alkalinity, and carbon isotope profiles in shallow sediment above marine gas hydrate systems: a numerical modeling perspective,” Journal of Geophysical Research, vol. 116, no. B9, 2011.

M. Luo, L. Chen, S. Wang, W. Yan, H. Wang, and D. Chen, “Pockmark activity inferred from pore water geochemistry in shallow sediments of the pockmark field in southwestern Xisha Uplift, northwestern South China Sea,” Marine and Petroleum Geology, vol. 48, pp. 247–259, 2013.

M. J. Whiticar, “Carbon and hydrogen isotope systematics of bacterial formation and oxidation of methane,” Chemical Geology, vol. 161, no. 1-3, pp. 291–314, 1999.

W. M. Sackett, “Carbon and hydrogen isotope effects during the thermocatalytic production of hydrocarbons in laboratory simulation experiments,” Geochimica et Cosmochimica Acta, vol. 42, no. 6, pp. 571–580, 1978.

W. S. Borowski, N. Cagatay, Y. Ternois, and C. K. Paull, “Data report: carbon isotopic composition of dissolved CO₂, CO₂ gas, and methane, Blake-Bahama Ridge and northeast Bermuda Rise, ODP Leg 172,” in Proceedings of the Ocean Drilling Program, 172 Scientific Results, L. D. Keigwin, D. Rio, G. D. Acton, and E. Arnold, Eds., pp. 1–16, ODP Publications, 2000.

Y. Chen, W. Ussler III, H. Haflidason et al., “Sources of methane inferred from pore-water δ¹³C of dissolved inorganic carbon in pockmark G11, offshore Mid-Norway,” Chemical Geology, vol. 275, no. 3-4, pp. 127–138, 2010.

Y. Hu, M. Luo, L. Chen et al., “Methane source linked to gas hydrate system at hydrate drilling areas of the South China Sea: porewater geochemistry and numerical model constraints,” Journal of Asian Earth Sciences, vol. 168, pp. 87–95, 2018.

F. Chen, J. Chen, H. Jin, and H. Li, “Correlation of δ¹³C origin surface sediments with sinking particulate matter in South China Sea and implication for reconstructing paleo- environment,” Acta Sedimentology Sinica, vol. 30, pp. 340–345, 2012.

X. Hu and D. J. Burdige, “Enriched stable carbon isotopes in the pore waters of carbonate sediments dominated by seagrasses: evidence for coupled carbonate dissolution and reprecipitation,” Geochimica et Cosmochimica Acta, vol. 71, no. 1, pp. 129–144, 2007.

X. Hu, W. J. Cai, Y. Wang, S. Luo, and X. Guo, “Pore-water geochemistry of two contrasting brine-charged seep sites in the northern Gulf of Mexico continental slope,” Marine Chemistry, vol. 118, no. 3-4, pp. 99–107, 2010.

G. E. Claypool and K. A. Kvenvolden, “Methane and other hydrocarbon gases in marine sediment,” Annual Review of Earth and Planetary Sciences, vol. 11, no. 1, pp. 299–327, 1983.

Q. Liang, S. Yang, J. Liang, and J. Tao, “The submarine seepage features of Haima cold seeps in the western part of the northern slope of the South China Sea,” in The 3rd Seafloor Observation Symposium, Qingdao, October 2016.

H. Guan, D. Birgel, J. Peckmann et al., “Lipid biomarker patterns of authigenic carbonates reveal fluid composition and seepage intensity at Haima cold seeps, South China Sea,” Journal of Asian Earth Sciences, vol. 168, pp. 163–172, 2018.