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 SO42-, CH4, PO43-, Ca2+, Mg2+, and dissolved inorganic carbon (DIC) in CL44 and CL47 using a steady-state reaction-transport model, we estimated that the dissolved SO42- was predominantly consumed by anaerobic oxidation of methane (AOM) at rates of 74.3 mmol m−2 yr−1 in CL44 and 85.0 mmol m−2 yr−1 in CL47. The relatively high AOM rates were sustained by free gas dissolution rather than local methanogenesis. Based on the diffusive Ba2+ 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 δ13CDIC mass balance calculation, we estimated that the δ13C 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 [912]. 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, Ba2+ 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 Ba2+ into the porewater below the SMTZ due to unsaturation. Through this cycling, authigenic barium fronts would be stably developed just above the SMTZ [1517]. The content of authigenic barite depends on the upward diffusive Ba2+ 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 Ba2+ 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 [1620].

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 [2126]. 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 [2731]. 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 (SO42-), calcium (Ca2+), magnesium (Mg2+), barium (Ba2+), phosphate (PO43-), methane (CH4), 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 δ13CDIC 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 [2731, 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).

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.

PO43- 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 cm3) onboard at core CL44. The porosity and density were calculated assuming a density of the porewater of 1.0 g cm-3.

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, SO42-, Ca2+, and Mg2+ 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%. Ba2+ concentrations were measured by inductively coupled plasma mass spectrometry (ICP-MS, Finnigan Element II). Before measurement, samples were prepared by diluting in 2% HNO3 with 10 ppb of Rh as an internal standard. The analytical precisions were estimated to be <5% for Ba2+. For cores CL30 and CL47, SO42-, Ca2+, and Mg2+ 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% HNO3 with 10 ppb of Tb as an internal standard before analysis for Ba2+ using the ICP-MS (Thermo Fisher iCAPQ). The analytical precisions were estimated to be <5% for Ba2+.

For core CL44, DIC concentrations and δ13CDIC values were determined using a continuous flow mass spectrometer (Thermo Fisher Delta-Plus). 0.5 ml porewater was treated with pure H3PO4 in a glass vial at 25°C. The CO2 produced was stripped with He and transferred into the mass spectrometer through the measurement of the δ13C 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 H3PO4 in a glass vial at 25°C. The CO2 produced was stripped with He and transferred into the mass spectrometer through which the δ13C values were measured. The analytical precisions were better than 0.2‰ for δ13C 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 HClO4 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 (CaCO3) 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% () H2O2 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 Ba2+ 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 Ba2+ (mmol m-2 yr-1), is the porosity, is the diffusion coefficient for seawater (m2 s-1), is the diffusion coefficient for sediments (m2 s-1), is the concentration of barium (mmol l-1), 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 Ba2+ 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 Ba2+. The diffusive flux was calculated using equations (1) and (2). is the tortuosity- and temperature-corrected diffusion coefficient of Ba2+ in the sediment, calculated from the diffusion coefficient in free solution () of 4.64, 4.62, and  cm2 s-1 (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 SO42-, CH4, DIC, PO43-, Ca2+, and Mg2+. The model is modified from previous simulations of methane-rich sediments [4851], 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.

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-3) and at depth below the sediment surface (mmol cm-3) (, Table 2). Bubble irrigation is described by parameters (yr-1) 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-3 yr-1), was described using a pseudo-first-order kinetic expression of the departure from the local methane gas solubility concentration, (mmol cm-3), where (yr-1) 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 CH2O and POP denote particulate organic carbon and phosphate, respectively. The total rate of POM mineralization, (wt.% C yr-1), 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 CO2 reduction [54]. Their net reactions at steady state are balanced with equivalent amounts of CO2 and CH4 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, (cm3 mmol-1 yr-1), is tuned to the sulfate profiles within the SMTZ.

The loss of Ca2+ and Mg2+ 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 CO32- 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 CH4. CH4 concentration at the lower boundary was a tunable parameter constrained from the SO42- profile. The model was solved using the NDSolve object of MATHEMATICA V. 10.0. The steady-state simulations were run for 107 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 SO42- concentration showed kink-type features at all the three cores (Figures 24 and Table 3). At site CL30, SO42- 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, SO42- 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). Ca2+ and Mg2+ 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. Ca2+ and Mg2+ concentrations declined sharply below the kinks due to ongoing carbonate precipitation and reached minimum at the SMTZ (Figures 24 and Table 3). Concentrations of DIC and PO43- showed opposite trends to SO42-, being depleted within the upper layer and enriched below it with the maximum at the STMZ (Figures 24 and Table 3). Moreover, CH4 concentrations at the three cores sharply increased below the SMTZ. The scatter in the CH4 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 δ13CDIC values (minimum values: -46.4‰ for CL30, -41.0‰ for CL44, and -38.8‰ for CL47) (Figures 24 and Table 3).

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 Ba2+ 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-1 (Table S6) with averages of 461 mg kg-1 for CL30, 502 mg kg-1 for CL44, and 502 mg kg-1 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-1 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-1, 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.

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 24 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-2 yr-1 for CL30, 1.54 mmol m-2 yr-1 for CL44, and 1.61 mmol m-2 yr-1 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, Ba2+ fluxes, and porosity are shown in Figures S5&S6.

4.3. Reaction-Transport Modeling

The modeled profiles and reaction rates are shown in Figures 24 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 SO42-, DIC, Ca2+, Mg2+, and PO43- at sites CL44, CL47, and CL30 above the kink with obvious discrepancies between modeled and measured concentrations of CH4 due to aforementioned degassing during core recovery (Figures 24). At site CL30, the model failed to reproduce the concentration gradients of SO42-, DIC, Ca2+, Mg2+, and PO43- below the kink (~3.5 mbsf) which is likely caused by a transient condition that is not considered in the steady-state model.

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 CH4 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 CH4 concentration reproduces the observed SO42- 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-2 y-1 () to 140 mmol m-2 y-1 ().

The initial age of the organic matter was tuned until a good fit was obtained for the PO43-. 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-2 yr-1) than that at site CL30 (18.8 mmol m-2 yr-1) (Table 4). The rates of POC degradation through sulfate reduction (POCSR) were 7.6, 23.9, and 25.1 mmol m-2 yr-1 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-2 yr-1 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 CH4 with benthic CH4 fluxes of 0.49, 2.0, and 2.7 mmol m-2 yr-1 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, 5964]. 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 CaCO3 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 [1013, 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 SO42- concentration profile can be reproduced after a sudden increase in CH4 concentration lasting for ~85 yr. On the other hand, based on the calculated diffusive Ba2+ 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 [7477] 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, 1619, 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 δ13CDIC 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 Ca2+ and Mg2+ concentrations with depth (Figures 36).

In addition, the δ13CDIC values below the SMTZ become more positive than those at the SMTZ. The reversal in δ13CDIC below the SMTZ is caused by the generation of 13C-enriched DIC via local methanogenesis at the methanogenic zone [63, 78]. The 13C-enriched DIC would migrate into the SMTZ from the methanogenic zone and “dilute” the 12C 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 CO2 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 δ13C values of methane are generally distinct between these two types of methane. The δ13C values of microbial methane typically range from −50‰ to −110‰ [80], whereas those of thermogenic methane range from −30‰ to −50‰ [81].

Because δ13C values of headspace methane in the sediments are absent in the study area, porewater DIC content and δ13CDIC 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) [8284]. In a closed system, δ13C 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 δ13CSW is assumed to be 0‰. The δ13C value of sedimentary organic matter in sediment of the SCS (−20‰; [85]) is used for the δ13COM. The overall δ13C of DIC derived from methanogenesis is equal to the δ13COM 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 δ13Cex can be acquired from a regression of porewater vs. DIC (Figure 7). The regression commonly shows linear in seep-impacted sediments, thus providing definitive δ13CDIC supplied to porewater [82, 86, 87]. The δ13Cex values and the contribution fractions estimated from the model are listed in Table 5. The estimated δ13CEM 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.

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 SO42-, CH4, DIC, PO43-, Ca2+, and Mg2+ 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−2 yr−1 for CL44 and 85.0 mmol m−2 yr−1 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 δ13C 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.


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 CaCO3 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 Ba2+ 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)