# An Areal Assessment of Subseafloor Carbon Cycling in Cold Seeps and Hydrate-Bearing Areas in the Northern South China Sea

Authors: Yanping Zhang; Min Luo; Yu Hu; Hongbin Wang; Duofu Chen
Journal: Geofluids
 (2019)

---

## Abstract

Gas hydrates, acting as a dynamic methane reservoir, store methane in the form of a solid phase under high-pressure and low-temperature conditions and release methane through the sediment column into seawater when they are decomposed. The seepage of methane-rich fluid (i.e., cold hydrocarbon seeps) fuels the chemosynthetic biota-inhabited surface sediments and represents the major pathway to transfer carbon from sediments to the water column. Generally, the major biogeochemical reactions related to carbon cycling in the anoxic marine sediments include organic matter degradation via sulfate reduction (OSR), anaerobic oxidation of methane (AOM), methanogenesis (ME), and carbonate precipitation (CP). In order to better understand the carbon turnover in the cold seeps and gas hydrate-bearing areas of the northern South China Sea (SCS), we collected geochemical data of 358 cores from published literatures and retrieved 37 cores and corresponding pore water samples from three areas of interest (i.e., Xisha, Dongsha, and Shenhu areas). Reaction-transport simulations indicate that the rates of organic matter degradation and carbonate precipitation are comparable in the three areas, while the rates of AOM vary over several orders of magnitude (AOM: 8.3-37.5 mmol·m-2·yr-1 in Dongsha, AOM: 12.4-170.6 mmol·m-2·yr-1 in Xisha, and AOM: 9.4-30.5 mmol·m-2·yr-1 in Shenhu). Both the arithmetical mean and interpolation mean of the biogeochemical processes were calculated in each area. Averaging these two mean values suggested that the rates of organic matter degradation in Dongsha (25.7 mmol·m-2·yr-1) and Xisha (25.1 mmol·m-2·yr-1) are higher than that in Shenhu (12 mmol·m-2·yr-1) and the AOM rate in Xisha (135.2 mmol·m-2·yr-1) is greater than those in Dongsha (27.8 mmol·m-2·yr-1) and Shenhu (17.5 mmol·m-2·yr-1). In addition, the rate of carbonate precipitation (32.3 mmol·m-2·yr-1) in Xisha is far higher than those of the other two regions (5.3 mmol·m-2·yr-1 in Dongsha, 5.8 mmol·m-2·yr-1 in Shenhu) due to intense AOM sustained by gas dissolution. In comparison with other cold seeps around the world, the biogeochemical rates in the northern SCS are generally lower than those in active continental margins and special environments (e.g., the Black sea) but are comparable with those in passive continental margins. Collectively, ~2.8 Gmol organic matter was buried and at least ~0.82 Gmol dissolved organic and inorganic carbon was diffused out of sediments annually. This may, to some extent, have an impact on the long-term deep ocean carbon cycle in the northern SCS.


---

## Body

## 1. Introduction

Marine sediments are the Earth’s largest methane reservoir, in which methane is dominantly preserved as gas hydrates that are commonly distributed along continental margins [1, 2]. It is estimated that at least 600 Gt methane is stored as hydrates in the continental margin sediments characterized by relatively high-pressure and low-temperature conditions [3]. However, as a dynamic methane reservoir, gas hydrates tend to decompose when the equilibrium state breaks, potentially leading to ocean acidification, submarine slope instability, and global climate change [2–4]. In addition, in light of the resource potential and environmental effect, gas hydrate has drawn increasing attentions worldwide since its first discovery in 1970s.Gas hydrates are usually closely coupled with cold seeps at continental slopes. On the continental margins, methane produced through either microbial methanogenesis or pyrolysis of organic matter can precipitate as hydrates when the dissolved methane concentration exceeds methane hydrate solubility within the gas hydrate stability zone [2, 5–7]. The gas hydrates will decompose and release large amounts of methane into the water column upon the pressure and/or temperature changes induced by global warming and sea level changes. The migration of methane-rich fluid towards the seafloor forms the cold seeps and associated seafloor expressions (e.g., mound or crater-like shape) [8]. Globally, the vast majority of methane in the sediments is consumed by the anaerobic oxidation of methane (AOM) mediated by anaerobic methane-oxidizing archaea and sulfate-reducing bacteria [9, 10]. AOM thus represents a microbial filter largely preventing dissolved methane from escaping the sediments. In the intensive seepage areas, a large amount of methane is usually transported in a gas phase as bubble, which cannot be fully consumed by AOM. The main product of AOM, dissolved inorganic carbon (DIC), can either precipitate as authigenic carbonate or migrate towards the seafloor. Other DIC sources in the sediments include organoclastic sulfate reduction and methanogenesis, which also makes important contribution to the DIC reservoir.The South China Sea (SCS), as one of the major marginal seas around the west Pacific Ocean, is well known for its extensive distribution of high-saturation gas hydrate [11–14]. Since the 1990’s, high-resolution seismic surveys have been carried out to delineate the distribution of bottom-simulating reflectors (BSR) and determine the prospecting area of gas hydrates. Through several hydrate survey expeditions, six promising hydrate-bearing areas including Dongsha, SW Taiwan, Xisha, Qiongdongnan, Shenhu, and Beikang have been confirmed to date [13–20]. Besides, more than 40 cold seepage sites have been indicated by the occurrence of 13C-depleted authigenic carbonate and seep-associated fauna and anomalous pore water and sediment geochemistry influenced by fluid seepage [6, 21–25]. Several pioneering studies have targeted at quantifying the rate of biogeochemical reactions and carbon fluxes at the seafloor using the reaction-transport model, in the specific seep sites of the northern SCS [16, 24, 26, 27]. Nevertheless, our quantitative understanding of the subsurface carbon cycle on a regional scale still remains limited. Therefore, the area-based carbon cycle calls for sufficient geochemical data combined with the numerical modeling approach to understand the relationship among different carbon reservoirs and their potential impact on the seawater carbonate system. In this study, 37 sediment cores were collected from the Shenhu area, Dongsha area, and Xisha area. The 1-D reaction-transport model was subsequently applied to quantify site-specific rates of biogeochemical processes and fluxes of DIC and methane. In combination with 358 published cores from the areas of interest, spatial interpolation was then used to explore the regional distribution of biogeochemical rates in the three regions. Finally, the areal assessments of surface carbon cycling in the northern SCS were provided to reveal the relationship among different carbon inventories.


## 2. Study Area

The SCS formed in the late Jurassic-early Cretaceous is located at the confluence of three plate collision, including the Pacific plate, Eurasia plate, and India-Australia plate [28]. The northern margin is a typical passive continental setting with a broad shelf, which is bounded to the west by the Indochina peninsula and to the east by a chain of island arcs [28]. The study area includes three areas in the northern slope of the SCS (Xisha, Shenhu, and Dongsha), with the water depth ranging from ~700 to ~2000 meters (see Figure 1). All these regions unanimously have experienced two tectonic stages, i.e., the syn-rift stage during the Paleocene and Oligocene and post-rift stage [29]. The infill of the northern continental marginal basins evolved simultaneously from alluvial and lacustrine deposits to neritic deposits, followed by progradational packages of slope sediments [30]. The lacustrine sediments deposited in the depressions acted as important source rocks for hydrocarbon generation. The potential abundant gas hydrate reservoirs have been indicated by the widespread occurrence of bottom-simulating reflectors (BSRs) [12, 31] (see Figure 1). The gas hydrates, retrieved by GMGS1&3 hydrate-drilling expeditions conducted by the Guangzhou Marine Geology Survey at the water depth of ~2000 meters in the Shenhu area, are characterized by a high-saturation (up to 40%) and mixed gas source [12]. During the GMGS2 hydrate-drilling expedition in Dongsha, various gas hydrate morphologies including massive, fracture-filling, and disseminated hydrates were found in the collected cores [31]. In addition, active seepage sites were discovered during remotely operated vehicle (ROV) deployments in the Dongsha and Xisha areas [21, 32]. These active seepage sites are characterized by the existence of gas ebullition in the water column; living chemosynthesis-based communities, such as mussels, tubeworms, and clams; and abundant authigenic carbonate rocks in the form of crusts, nodules, and tubular concretion.Figure 1
The bathematric map of the northern South China Sea. Rectangle boxes denote study areas. Note that Xisha herein includes the Qiongdongnan basin and Xisha Uplift, and Dongsha comprises a broad area of the eastern Zhujiangkou basin and SW Taiwan basin. The blue stars denote active seepage sites. The red triangle represents the Jiulong methane reef covering an area of ~400 km2.


## 3. Sampling and Methods

### 3.1. Sampling

In total, the 37 sediment cores used in this study were collected from five expeditions by using different research vessels. Four cores were retrieved by gravity corers from the Dongsha area during the Haiyang-4 cruise conducted by the Guangzhou Marine Geological Survey in October 2013. The sediments consisted of homogeneous green-gray clay silt, embedded by several small authigenic carbonate nodules. The Dongsha area was visited twice by Shiyan-3 in April 2014 and by Haiyang-4 in May 2015. The sediments were sampled by piston-gravity corers and box corers. In the Xisha area, five piston cores were retrieved by Haiyang-4 in April 2015. The sediments change from brown silt clay to gray calcareous clay. Twenty-four cores in Shenhu were collected by piston corer using Haiyang-4 in May 2015 and September 2016. The lithology in the Shenhu area was similar to that in the Xisha area, consisting of brown to gray calcareous silt. Detailed information with regard to core lengths, water depths, and sampling gears were shown in TableS1. All the retrieved cores were immediately brought to the onboard laboratory for porewater extraction. Rhizon samplers were used to collect porewater at 20 cm or 40 cm intervals. Porewater aliquots were stored at 4°C with ultrapure concentrated HNO3 and saturated HgCl2 solution for anion and cation analyses, respectively.


### 3.2. Analytical Methods

Concentrations of total alkalinity (TA) and phosphate (PO43-) were measured onboard immediately after pore water collection. The TA concentrations were determined by direct titration with 0.006 M HCl using Bruevich’s method [33]. The analysis was calibrated using standard seawater (IAPSO) with a precision of better than 2%. Phosphate concentrations were determined using a HITACHI U-5100 spectrophotometer with the analytical precision of better than 5%. Concentrations of sulfate (SO42-), chloride (Cl-), and calcium (Ca2+) were determined onshore by a Dionex ICS-5000+ ion chromatography at the South China Sea Institute of Oceanology, Chinese Academy of Sciences (CAS), with the analytical precision better than 2%. Dissolved inorganic carbon (DIC) concentrations were measured by an IsoPrime 100 continuous flow isotope ration mass spectrometer (CF-IRMS) at the State Key Laboratory of Isotope Geochemistry, Guangzhou Institute of Geochemistry, CAS, with the precision of better than 2%. Detailed analytical processes can be found in the study of Hu et al. [34].


### 3.3. Reaction-Transport Modeling

A one-dimensional, steady-state, reaction-transport model developed from previous approaches was applied to simulate one as a solid phase (POC) and five as dissolved species (SO42-, CH4, DIC, Ca2+, and PO43-). Two partial differential equations were used to reproduce the depth concentration profiles of solid and dissolved species [35–37].
(1)ϕ∂Ca∂t=∂ϕ·Ds·∂Ca/∂x∂x−∂ϕ·vp·Ca∂x+ϕ·∑R+ϕ·RBui,(2)1−ϕ∂POC∂t=−∂1−ϕ·vs·POC∂x+1−ϕ·∑R,where t (year) is time, x (cm) is depth below the seafloor, ϕ (dimensionless) is porosity, Ds (cm2·yr-1) is solute-molecular diffusion coefficient corrected by tortuosity, Ca (mol/L) is the concentrations of solutes in pore water, POC (wt.%) is the content of POC in dry sediment, vp (cm·yr-1) is the burial velocity of porewater due to steady-state compaction, vs (cm·yr-1) is burial velocity of solid, ∑R is the sum of rates of all chemical reactions considered in the model, and RBui is the mixing rate of bottom water and porewater due to bubble irrigation.Sediment porosity decreases with depth assuming steady-state compaction.vp and vs are variables changing with porosity under assumption of steady-state compaction. Equations can be described as follows:
(3)ϕ=ϕf+ϕ0−ϕf·e−px,where ϕf (dimensionless) is the porosity at depth and ϕ0 (dimensionless) is the porosity at the sediment-water interface. p (cm-1) denotes the attenuation coefficient of porosity.In the absence of externally imposed fluid advection at the seafloor, the velocity of porewater and solids is directed downward under steady-state compaction relative to the seafloor:(4)vp=ϕf·wϕ,vs=1−ϕf·w1−ϕ,where w (cm·yr-1) is the sedimentation velocity.Calibration of diffusive coefficient with tortuosity was after the equation of Boudreau [36].
(5)Ds=Dm1−lnϕ2,where Dm (cm2·yr-1) is molecular diffusion coefficient in the in situ temperature, salinity, and pressure calculated according to Li and Gregory [38].In sites 2015XS-R2 and 2015XS-50, gas ebullition during the ROV deployment was observed [32]. The porewater mixing with bottom water induced by rising gas bubbles can be described as a nonlocal transport similar to bioirrigation [39]:
(6)RBui=α0·expLirr−x/α11+expLirr−x/α1·C0−Cl,where α0 (yr-1) is the coefficient of irrigation intensity, Lirr (cm) is the depth of bubble irrigation, α1 (cm) is the parameter determining how quickly bubble irrigation is attenuated to zero at an approximate depth of Lirr, C0 is solute concentration at the SWI, and Cl is concentration at any depth within the irrigation zone.As for biogeochemical reactions, organic matter decomposition via sulfate (OSR), methanogenesis (ME), anaerobic oxidation of methane (AOM), and authigenic carbonate precipitation (CP), as well as methane gas dissolution in a shallow hydrate-bearing site, are included in the model. Aerobic respiration, denitrification, manganese reduction, and iron reduction via organic matter remineralization were not taken into account in the model since these reactions were only restricted to the uppermost sediments (about 10-20 cm below the seafloor). Below the oxic and suboxic zones, sulfate serves as the main terminal electron acceptor for oxidizing organic matter [40–43]. Organic matter degradation via sulfate reduction was expressed as follows:
(7)2CH2OPO43−rp+SO42−→2HCO3−+H2S+2rpPO43−,where rp is the Redfield ratio of organic phosphorus to organic carbon.Below the sulfate reduction zone, the organic matter was degraded via methanogenesis:(8)2CH2OPO43−rp→CO2+CH4+2rpPO43−.The rates of sulfate reduction and methanogenesis depend on the total rate of POC degradation according to Middelburg [44] and Wallmann et al. [37]:
(9)RPOC=KcDIC+CH4+Kc·0.16·a0+xw−0.95·POC,where Kc (M) is an inhibition coefficient of POC degradation, DIC and CH4 (mM) are concentrations of DIC and CH4, respectively, and a0 (yr) is the initial age of organic matter, which was constrained using measured PO43- concentrations. The rates of OSR and ME are thus expressed using the following equations:
(10)ROSR=0.5·ρs·1−ϕ·106MWc·ϕ·SO42−SO42−+KSO42−·RPOC,RME=0.5·ρs·1−ϕ·106MWc·ϕ·KSO42−SO42−+KSO42−·RPOC,where ρs (g·cm-3) is the density of dry sediments, MWc (g·mol-1) is the molecular weight of carbon, and KSO42− is the Michaelis-Menten constant for the inhibition of sulfate concentration.The AOM-coupled sulfate reduction mediated by a syntrophic consortium of methanotrophic archaea and sulfate-reducing bacteria is an important sink for methane in anoxic marine sediments [9, 45].
(11)CH4+SO42−→HCO3−+HS−+H2O.The rate ofAOM was calculated using bimolecular kinetics:
(12)RAOM=kAOM·SO42−·CH4,where kAOM (cm3·mol-1·yr-1) is an apparent second-order rate constant.The DIC produced by AOM and organic matter degradation can be consumed due to precipitation of authigenic carbonate [35, 46, 47]:
(13)Ca2++2HCO3−→CaCO3+CO2+H2O.The rate of authigenic carbonate precipitation was simulated using the thermodynamic solubility constant as defined by Millero [48]:
(14)RCP=kCP·Ca2+·CO32−KSP−1,where kCP (mol·cm-3·yr-1) is the kinetic constant and KSP (mol2·l2) is the thermodynamic equilibrium constant. A typical porewater pH value of 7.6 was used to calculate CO32- from modeled DIC concentration [49]. CaCO3 is not simulated explicitly in the model.At sites 2015XS-R2, 2015XS-50 where gas ebullition was observed, dissolution of rising gas bubble occurs if the porewater is undersaturated with respect to in situ methane gas solubility:(15)CH4g→CH4aq,whenCH4aq≤LMB,where LMB (mol/L) is a site-specific constant dependent on in situ salinity, pressure, and temperature using the algorithm in the work of Duan et al. [7, 50]. The rate of gas dissolution was described using a first-order kinetic expression of the departure from the local methane gas solubility concentration:
(16)RMB=kMB·LMB−CH4,where kMB (yr-1) is a first-order rate constant.Gas hydrate precipitation occurred if the dissolved methane concentration exceeded the solubility of gas hydrate (LGH):
(17)RGH=kGH·CH4LGH−1,where kGH (mol·cm-3·yr-1) is a fitting parameter constrained by the depth of the first occurrence of gas hydrate.Upper and lower boundary conditions for all species were imposed as fixed concentrations with the exception of the intensive seepage site (2015XS-R2) whose lower boundary was imposed as a zero-concentration gradient. The continuous differential equations in equations (1) and (2) were solved using finite differences and the method-of-lines over an uneven grid with a higher spatial resolution at the surface and lower resolution towards the bottom. The model was solved using the NDSlove object of MATHEMATICA version 9.0. All simulations were run for long enough to achieve the steady state with a mass conservation of >99%. Parameter information and reaction terms of all species used in the model can be found in Table S2, Table S3, and Table S4.


### 3.4. Spatial Interpolation

Site-specific fluxes and rates were then used to explore area-based distribution in spatial coordinates by using the ordinary Kriging method. It is more accurate and realistic in this method to predict variables in the confined zone combining spatial distance with coordinate. According to Matheron [51], the mathematical principle of this method can be described as
(18)Z′x0=∑i=1nλiZxi,where Zxi is a value which is situated at any place of the study area, x0 is the predicted site, and λi is the weight factor, of which the sum is always equal to unity.Spatial interpolation was implemented by ARCMAP version 10.5. Interpolation extent was derived by the point feature. The cell size of the output raster which was automatically created was 1/250 of the width and height of the extent.


## 3.1. Sampling

In total, the 37 sediment cores used in this study were collected from five expeditions by using different research vessels. Four cores were retrieved by gravity corers from the Dongsha area during the Haiyang-4 cruise conducted by the Guangzhou Marine Geological Survey in October 2013. The sediments consisted of homogeneous green-gray clay silt, embedded by several small authigenic carbonate nodules. The Dongsha area was visited twice by Shiyan-3 in April 2014 and by Haiyang-4 in May 2015. The sediments were sampled by piston-gravity corers and box corers. In the Xisha area, five piston cores were retrieved by Haiyang-4 in April 2015. The sediments change from brown silt clay to gray calcareous clay. Twenty-four cores in Shenhu were collected by piston corer using Haiyang-4 in May 2015 and September 2016. The lithology in the Shenhu area was similar to that in the Xisha area, consisting of brown to gray calcareous silt. Detailed information with regard to core lengths, water depths, and sampling gears were shown in TableS1. All the retrieved cores were immediately brought to the onboard laboratory for porewater extraction. Rhizon samplers were used to collect porewater at 20 cm or 40 cm intervals. Porewater aliquots were stored at 4°C with ultrapure concentrated HNO3 and saturated HgCl2 solution for anion and cation analyses, respectively.


## 3.2. Analytical Methods

Concentrations of total alkalinity (TA) and phosphate (PO43-) were measured onboard immediately after pore water collection. The TA concentrations were determined by direct titration with 0.006 M HCl using Bruevich’s method [33]. The analysis was calibrated using standard seawater (IAPSO) with a precision of better than 2%. Phosphate concentrations were determined using a HITACHI U-5100 spectrophotometer with the analytical precision of better than 5%. Concentrations of sulfate (SO42-), chloride (Cl-), and calcium (Ca2+) were determined onshore by a Dionex ICS-5000+ ion chromatography at the South China Sea Institute of Oceanology, Chinese Academy of Sciences (CAS), with the analytical precision better than 2%. Dissolved inorganic carbon (DIC) concentrations were measured by an IsoPrime 100 continuous flow isotope ration mass spectrometer (CF-IRMS) at the State Key Laboratory of Isotope Geochemistry, Guangzhou Institute of Geochemistry, CAS, with the precision of better than 2%. Detailed analytical processes can be found in the study of Hu et al. [34].


## 3.3. Reaction-Transport Modeling

A one-dimensional, steady-state, reaction-transport model developed from previous approaches was applied to simulate one as a solid phase (POC) and five as dissolved species (SO42-, CH4, DIC, Ca2+, and PO43-). Two partial differential equations were used to reproduce the depth concentration profiles of solid and dissolved species [35–37].
(1)ϕ∂Ca∂t=∂ϕ·Ds·∂Ca/∂x∂x−∂ϕ·vp·Ca∂x+ϕ·∑R+ϕ·RBui,(2)1−ϕ∂POC∂t=−∂1−ϕ·vs·POC∂x+1−ϕ·∑R,where t (year) is time, x (cm) is depth below the seafloor, ϕ (dimensionless) is porosity, Ds (cm2·yr-1) is solute-molecular diffusion coefficient corrected by tortuosity, Ca (mol/L) is the concentrations of solutes in pore water, POC (wt.%) is the content of POC in dry sediment, vp (cm·yr-1) is the burial velocity of porewater due to steady-state compaction, vs (cm·yr-1) is burial velocity of solid, ∑R is the sum of rates of all chemical reactions considered in the model, and RBui is the mixing rate of bottom water and porewater due to bubble irrigation.Sediment porosity decreases with depth assuming steady-state compaction.vp and vs are variables changing with porosity under assumption of steady-state compaction. Equations can be described as follows:
(3)ϕ=ϕf+ϕ0−ϕf·e−px,where ϕf (dimensionless) is the porosity at depth and ϕ0 (dimensionless) is the porosity at the sediment-water interface. p (cm-1) denotes the attenuation coefficient of porosity.In the absence of externally imposed fluid advection at the seafloor, the velocity of porewater and solids is directed downward under steady-state compaction relative to the seafloor:(4)vp=ϕf·wϕ,vs=1−ϕf·w1−ϕ,where w (cm·yr-1) is the sedimentation velocity.Calibration of diffusive coefficient with tortuosity was after the equation of Boudreau [36].
(5)Ds=Dm1−lnϕ2,where Dm (cm2·yr-1) is molecular diffusion coefficient in the in situ temperature, salinity, and pressure calculated according to Li and Gregory [38].In sites 2015XS-R2 and 2015XS-50, gas ebullition during the ROV deployment was observed [32]. The porewater mixing with bottom water induced by rising gas bubbles can be described as a nonlocal transport similar to bioirrigation [39]:
(6)RBui=α0·expLirr−x/α11+expLirr−x/α1·C0−Cl,where α0 (yr-1) is the coefficient of irrigation intensity, Lirr (cm) is the depth of bubble irrigation, α1 (cm) is the parameter determining how quickly bubble irrigation is attenuated to zero at an approximate depth of Lirr, C0 is solute concentration at the SWI, and Cl is concentration at any depth within the irrigation zone.As for biogeochemical reactions, organic matter decomposition via sulfate (OSR), methanogenesis (ME), anaerobic oxidation of methane (AOM), and authigenic carbonate precipitation (CP), as well as methane gas dissolution in a shallow hydrate-bearing site, are included in the model. Aerobic respiration, denitrification, manganese reduction, and iron reduction via organic matter remineralization were not taken into account in the model since these reactions were only restricted to the uppermost sediments (about 10-20 cm below the seafloor). Below the oxic and suboxic zones, sulfate serves as the main terminal electron acceptor for oxidizing organic matter [40–43]. Organic matter degradation via sulfate reduction was expressed as follows:
(7)2CH2OPO43−rp+SO42−→2HCO3−+H2S+2rpPO43−,where rp is the Redfield ratio of organic phosphorus to organic carbon.Below the sulfate reduction zone, the organic matter was degraded via methanogenesis:(8)2CH2OPO43−rp→CO2+CH4+2rpPO43−.The rates of sulfate reduction and methanogenesis depend on the total rate of POC degradation according to Middelburg [44] and Wallmann et al. [37]:
(9)RPOC=KcDIC+CH4+Kc·0.16·a0+xw−0.95·POC,where Kc (M) is an inhibition coefficient of POC degradation, DIC and CH4 (mM) are concentrations of DIC and CH4, respectively, and a0 (yr) is the initial age of organic matter, which was constrained using measured PO43- concentrations. The rates of OSR and ME are thus expressed using the following equations:
(10)ROSR=0.5·ρs·1−ϕ·106MWc·ϕ·SO42−SO42−+KSO42−·RPOC,RME=0.5·ρs·1−ϕ·106MWc·ϕ·KSO42−SO42−+KSO42−·RPOC,where ρs (g·cm-3) is the density of dry sediments, MWc (g·mol-1) is the molecular weight of carbon, and KSO42− is the Michaelis-Menten constant for the inhibition of sulfate concentration.The AOM-coupled sulfate reduction mediated by a syntrophic consortium of methanotrophic archaea and sulfate-reducing bacteria is an important sink for methane in anoxic marine sediments [9, 45].
(11)CH4+SO42−→HCO3−+HS−+H2O.The rate ofAOM was calculated using bimolecular kinetics:
(12)RAOM=kAOM·SO42−·CH4,where kAOM (cm3·mol-1·yr-1) is an apparent second-order rate constant.The DIC produced by AOM and organic matter degradation can be consumed due to precipitation of authigenic carbonate [35, 46, 47]:
(13)Ca2++2HCO3−→CaCO3+CO2+H2O.The rate of authigenic carbonate precipitation was simulated using the thermodynamic solubility constant as defined by Millero [48]:
(14)RCP=kCP·Ca2+·CO32−KSP−1,where kCP (mol·cm-3·yr-1) is the kinetic constant and KSP (mol2·l2) is the thermodynamic equilibrium constant. A typical porewater pH value of 7.6 was used to calculate CO32- from modeled DIC concentration [49]. CaCO3 is not simulated explicitly in the model.At sites 2015XS-R2, 2015XS-50 where gas ebullition was observed, dissolution of rising gas bubble occurs if the porewater is undersaturated with respect to in situ methane gas solubility:(15)CH4g→CH4aq,whenCH4aq≤LMB,where LMB (mol/L) is a site-specific constant dependent on in situ salinity, pressure, and temperature using the algorithm in the work of Duan et al. [7, 50]. The rate of gas dissolution was described using a first-order kinetic expression of the departure from the local methane gas solubility concentration:
(16)RMB=kMB·LMB−CH4,where kMB (yr-1) is a first-order rate constant.Gas hydrate precipitation occurred if the dissolved methane concentration exceeded the solubility of gas hydrate (LGH):
(17)RGH=kGH·CH4LGH−1,where kGH (mol·cm-3·yr-1) is a fitting parameter constrained by the depth of the first occurrence of gas hydrate.Upper and lower boundary conditions for all species were imposed as fixed concentrations with the exception of the intensive seepage site (2015XS-R2) whose lower boundary was imposed as a zero-concentration gradient. The continuous differential equations in equations (1) and (2) were solved using finite differences and the method-of-lines over an uneven grid with a higher spatial resolution at the surface and lower resolution towards the bottom. The model was solved using the NDSlove object of MATHEMATICA version 9.0. All simulations were run for long enough to achieve the steady state with a mass conservation of >99%. Parameter information and reaction terms of all species used in the model can be found in Table S2, Table S3, and Table S4.


## 3.4. Spatial Interpolation

Site-specific fluxes and rates were then used to explore area-based distribution in spatial coordinates by using the ordinary Kriging method. It is more accurate and realistic in this method to predict variables in the confined zone combining spatial distance with coordinate. According to Matheron [51], the mathematical principle of this method can be described as
(18)Z′x0=∑i=1nλiZxi,where Zxi is a value which is situated at any place of the study area, x0 is the predicted site, and λi is the weight factor, of which the sum is always equal to unity.Spatial interpolation was implemented by ARCMAP version 10.5. Interpolation extent was derived by the point feature. The cell size of the output raster which was automatically created was 1/250 of the width and height of the extent.


## 4. Results

### 4.1. Porewater Geochemistry

We analysed 37 sediment cores, including 10 in Dongsha, 3 in Xisha, and 24 in Shenhu, and all the geochemical data are shown in FigureS1.The sulfate concentration profiles in Dongsha generally show quasilinear depletion with depth. Especially in 2013DS-F, 2014DS-G4, and 2015DS17, sulfate concentrations show sharper decline towards the sulfate-methane transition zone (SMTZ) compared with those in other sites. Calcium concentration profiles in Dongsha display linear decline with depth. In contrast, DIC concentration profiles in Dongsha generally mirror the calcium concentration profiles.The sulfate concentrations in Shenhu basically quasi-linearly decline with depth but their gradients differ among sites. Similar to the calcium concentration profiles in Dongsha, calcium concentration profiles in Shenhu display linear decrease with depth and DIC concentration profiles show quasi-linear increase with depth.The sulfate concentration profiles in Xisha do not follow the general trend of quasi-linear decrease. Instead, they show kink-type features with inflection points being located at 2 mbsf to 7 mbsf. The SO42- concentrations in 2015XS-50 and 2015XS-R2 display near-seawater values in the upper 3-5 meters and then decline sharply towards the SMTZ. The seawater intrusion feature is also reflected by calcium concentration profiles, which is attributed to the irrigation of porewater caused by gas bubbling as demonstrated by previous studies [26, 52, 53]. At core 2015XS-44, the sulfate concentration profile does not show the feature of bubble irrigation but the calcium concentration profile still exhibits a kink at ~5 mbsf. The reason is not well understood at this moment, but it might be related to carbonate dissolution induced by organoclastic aerobic oxidation. At cores 2015XS-50 and 2015XS-R2, DIC concentration profiles also show the bubble irrigation feature and generally mirror the trend of sulfate concentration profiles.


### 4.2. Reaction-Transport Modeling

The simulation results for the three areas mostly reproduced the measured profiles. Model parameters used to constrain the model curve are listed in TableS2 and Table S3. In this model, phosphate concentration profiles were used to constrain the initial age of POC. Due to the lack of phosphate concentration profiles in the cores of Dongsha, initial ages in Dongsha (5000 yr) were adopted from Hu et al. [54]. For the same reason, the initial ages of 2015XS-44 and 2015XS-50 in Xisha were assumed to be the same as that of 2015XS-R2 derived from a measured phosphate concentration profile. Likewise, in the Shenhu area, the average initial age of TOC for the cores with phosphate data (200 ka) was taken to model those cores that lack phosphate sites. The modeled POC degradation rates are 11.2-24.5 mmol·m-2·yr-1 in Dongsha, 2.9-7.5 mmol·m-2·yr-1 in Xisha, and 7.2-30.7 mmol·m-2·yr-1 in Shenhu.The depth-integrated rates and fluxes for the individual core are presented in TableS1. The rates of AOM in Dongsha range from 8.3 to 37.5 mmol·m-2·yr-1, which are similar to AOM rates in Shenhu (9.4 to 30.5 mmol·m-2·yr-1). In contrast, AOM rates in Xisha (12.4-170.6 mmol·m-2·yr-1) are greatly higher than those in Dongsha and Shenhu, because AOM was primarily supplied by methane gas dissolution. The methane fluxes at the top of the simulated sediment column are 0.002-0.2 mmol·m-2·yr-1 in Dongsha, 0.002-0.3 mmol·m-2·yr-1 in Shenhu, and 0.01-399.2 mmol·m-2·yr-1 in Xisha. DIC fluxes in all three areas are approximately three to four orders of magnitude higher than methane fluxes with 13.1-26.1 mmol·m-2·yr-1 in Dongsha, 10.1-31.7 mmol·m-2·yr-1 in Shenhu, and 1-155 mmol·m-2·yr-1 in Xisha. Rates of carbonate precipitation in Dongsha (4.3-7.3 mmol·m-2·yr-1), Shenhu (3.6-9 mmol·m-2·yr-1), and Xisha (4.1-12.4 mmol·m-2·yr-1) are comparable.


### 4.3. Spatial Interpolation

Interpolation was utilized to explore the spatial distribution of the biogeochemical rates of 395 cores, of which 37 cores are new in this study (see TableS1) and 358 cores were collected from literatures (see Table S5). As shown in Figure 2, rates of POC degradation via sulfate reduction and methanogenesis in Dongsha generally show high values in the east and gradually decrease towards the west. The rates of authigenic carbonate precipitation are mostly in agreement with those of AOM with higher values towards the east even though within a smaller region. The POC degradation rates in Xisha decrease from the south to the north. In contrast to that in the Dongsha area, the distribution of AOM rates in the Xisha are is inconsistent with that of carbonate precipitation rates. This may be due to the lack of coordinates in some cores with high carbonate precipitation rates which are not incorporated in the north. Unlike the POC decomposition rate distribution in Dongsha and Xisha, two spots of intensified POC decomposition are located to the northeast and to the southwest in Shenhu. Rates of AOM and carbonate precipitation are generally low in the whole area.Figure 2
Spatial distribution of rates of organoclastic sulfate reduction (OSR), anaerobic oxidation of methane (AOM), methanogenesis (ME), and authigenic carbonate precipitation (CP) in the three areas. The black dots represent spatial distribution of sampling sites. Different colors denote different ranges of biogeochemical rates. In comparison to the interpolation rates of OSR, ME, and AOM in Dongsha, a smaller area was confined in the interpolation rate of CP. Additionally, rates of OSR, ME, AOM, and CP have the same interpolation region and so do those in Xisha.


## 4.1. Porewater Geochemistry

We analysed 37 sediment cores, including 10 in Dongsha, 3 in Xisha, and 24 in Shenhu, and all the geochemical data are shown in FigureS1.The sulfate concentration profiles in Dongsha generally show quasilinear depletion with depth. Especially in 2013DS-F, 2014DS-G4, and 2015DS17, sulfate concentrations show sharper decline towards the sulfate-methane transition zone (SMTZ) compared with those in other sites. Calcium concentration profiles in Dongsha display linear decline with depth. In contrast, DIC concentration profiles in Dongsha generally mirror the calcium concentration profiles.The sulfate concentrations in Shenhu basically quasi-linearly decline with depth but their gradients differ among sites. Similar to the calcium concentration profiles in Dongsha, calcium concentration profiles in Shenhu display linear decrease with depth and DIC concentration profiles show quasi-linear increase with depth.The sulfate concentration profiles in Xisha do not follow the general trend of quasi-linear decrease. Instead, they show kink-type features with inflection points being located at 2 mbsf to 7 mbsf. The SO42- concentrations in 2015XS-50 and 2015XS-R2 display near-seawater values in the upper 3-5 meters and then decline sharply towards the SMTZ. The seawater intrusion feature is also reflected by calcium concentration profiles, which is attributed to the irrigation of porewater caused by gas bubbling as demonstrated by previous studies [26, 52, 53]. At core 2015XS-44, the sulfate concentration profile does not show the feature of bubble irrigation but the calcium concentration profile still exhibits a kink at ~5 mbsf. The reason is not well understood at this moment, but it might be related to carbonate dissolution induced by organoclastic aerobic oxidation. At cores 2015XS-50 and 2015XS-R2, DIC concentration profiles also show the bubble irrigation feature and generally mirror the trend of sulfate concentration profiles.


## 4.2. Reaction-Transport Modeling

The simulation results for the three areas mostly reproduced the measured profiles. Model parameters used to constrain the model curve are listed in TableS2 and Table S3. In this model, phosphate concentration profiles were used to constrain the initial age of POC. Due to the lack of phosphate concentration profiles in the cores of Dongsha, initial ages in Dongsha (5000 yr) were adopted from Hu et al. [54]. For the same reason, the initial ages of 2015XS-44 and 2015XS-50 in Xisha were assumed to be the same as that of 2015XS-R2 derived from a measured phosphate concentration profile. Likewise, in the Shenhu area, the average initial age of TOC for the cores with phosphate data (200 ka) was taken to model those cores that lack phosphate sites. The modeled POC degradation rates are 11.2-24.5 mmol·m-2·yr-1 in Dongsha, 2.9-7.5 mmol·m-2·yr-1 in Xisha, and 7.2-30.7 mmol·m-2·yr-1 in Shenhu.The depth-integrated rates and fluxes for the individual core are presented in TableS1. The rates of AOM in Dongsha range from 8.3 to 37.5 mmol·m-2·yr-1, which are similar to AOM rates in Shenhu (9.4 to 30.5 mmol·m-2·yr-1). In contrast, AOM rates in Xisha (12.4-170.6 mmol·m-2·yr-1) are greatly higher than those in Dongsha and Shenhu, because AOM was primarily supplied by methane gas dissolution. The methane fluxes at the top of the simulated sediment column are 0.002-0.2 mmol·m-2·yr-1 in Dongsha, 0.002-0.3 mmol·m-2·yr-1 in Shenhu, and 0.01-399.2 mmol·m-2·yr-1 in Xisha. DIC fluxes in all three areas are approximately three to four orders of magnitude higher than methane fluxes with 13.1-26.1 mmol·m-2·yr-1 in Dongsha, 10.1-31.7 mmol·m-2·yr-1 in Shenhu, and 1-155 mmol·m-2·yr-1 in Xisha. Rates of carbonate precipitation in Dongsha (4.3-7.3 mmol·m-2·yr-1), Shenhu (3.6-9 mmol·m-2·yr-1), and Xisha (4.1-12.4 mmol·m-2·yr-1) are comparable.


## 4.3. Spatial Interpolation

Interpolation was utilized to explore the spatial distribution of the biogeochemical rates of 395 cores, of which 37 cores are new in this study (see TableS1) and 358 cores were collected from literatures (see Table S5). As shown in Figure 2, rates of POC degradation via sulfate reduction and methanogenesis in Dongsha generally show high values in the east and gradually decrease towards the west. The rates of authigenic carbonate precipitation are mostly in agreement with those of AOM with higher values towards the east even though within a smaller region. The POC degradation rates in Xisha decrease from the south to the north. In contrast to that in the Dongsha area, the distribution of AOM rates in the Xisha are is inconsistent with that of carbonate precipitation rates. This may be due to the lack of coordinates in some cores with high carbonate precipitation rates which are not incorporated in the north. Unlike the POC decomposition rate distribution in Dongsha and Xisha, two spots of intensified POC decomposition are located to the northeast and to the southwest in Shenhu. Rates of AOM and carbonate precipitation are generally low in the whole area.Figure 2
Spatial distribution of rates of organoclastic sulfate reduction (OSR), anaerobic oxidation of methane (AOM), methanogenesis (ME), and authigenic carbonate precipitation (CP) in the three areas. The black dots represent spatial distribution of sampling sites. Different colors denote different ranges of biogeochemical rates. In comparison to the interpolation rates of OSR, ME, and AOM in Dongsha, a smaller area was confined in the interpolation rate of CP. Additionally, rates of OSR, ME, AOM, and CP have the same interpolation region and so do those in Xisha.


## 5. Discussion

### 5.1. Regional Biogeochemical Rates and Fluxes of Methane and DIC

By upscaling the rates and fluxes in the specific sites into a region of interest, raster means were obtained by ARCMAP based on the small subareas divided automatically, which represent areal means. To better constrain the areal mean values of individual biogeochemical reactions in each area, we averaged arithmetical mean and raster mean values (see Table1). Among the three areas, standard deviations of arithmetical means in the Xisha area were generally higher than 100, which indicate high heterogeneity of the biogeochemical rates. The arithmetical average of the rates of AOM in the Xisha area was much greater than those in the other two areas due to intensive methane bubbling. However, the raster mean of AOM rates in Xisha was slightly lower than those at the other two areas. The inconsistent result derived using two different methods is ascribed to the exclusion of intensive seepage sites due to unavailability of the sampling locations from the interpolation. Regardless of all the discrepancies between the two methods, a tentative comparison among the three regions indicates that the rates of AOM and carbonate precipitation in Xisha are overwhelmingly larger than those in the other two areas and that the rate of organic matter degradation in Xisha was slightly higher than those in Dongsha and Shenhu areas. The low rates of biogeochemical reactions in Shenhu are likely due to the absence of methane seepage and the occurrence of deep-seated hydrates [12]. The Dongsha area always shows the intermediate rates compared with the other two areas owing to the occurrence of sporadic methane seepage sites and gas hydrates in shallow sediments.Table 1
Biogeochemical rates and effluxes of methane and DIC in the three areas.
Location
CH4 efflux
DIC efflux
OSR
ME
AOM
CP
Cores
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
Dongsha
0.2
9.7
95
4.9
25.3
19.3
7
22.3
14.1
17.7
12
15.9
5.1
4.8
4
4.9
29.0
26.7
14
27.8
5.5
5.4
1
5.3
363
Xisha
0.1
144
366
72
36.0
246.2
136
141
27.6
14.6
20
21.1
2.6
1.5
1
2
9.4
261.1
385
135.2
27.3
37.4
28
32.3
8
Shenhu
0.0
0.1
1
0.0
17.1
19.8
5
18.4
5.8
9.8
2
7.8
1.4
2.8
1
2.1
16.8
18.2
5
17.5
5.4
6.3
1.6
5.8
24
North SCS
0.13
11.3
96
5.7
29.5
62.9
178
46.2
19.6
13.4
12
16.5
3.6
3.2
3
3.4
19.0
60.8
180
39.9
15.4
12
17
13.7
395
RM: raster mean; AM: arithmetical mean; SD: standard deviation.The Xisha area is characterized by high content of sedimentary organic matter (up to 2%) mainly sourced from marine algae [55, 56]. This may explain the relatively higher rate of organic matter degradation in Xisha compared to Dongsha and Shenhu. Applying the bimolecular kinetics to calculate the AOM rate, we assume that the adjustable rate constant (kAOM) represents the integrated effects of a number of factors including the microbial community structure and abundance, bioenergetics, and enzyme kinetics. Therefore, different settings may have distinct rate constants, ranging over 6 orders of magnitude [57]. Similar to the other passive continental margin environments, the rate constant (kAOM=1) in this study is generally lower than those in active continental margin environments [26, 37, 58–61]. In the cold seeps and hydrate-bearing area, methane expulsion is usually associated with some specific structures, e.g., mud diapir, fault, and gas chimney. These structures provide conduits for the upward migration of methane-rich fluid. According to the geophysical imaging in Dongsha, large amounts of mud diapir and gas chimney are well developed in subsurface sediments [62]. Gas bubble release was indicated by hydroacoustic flares, and indeed, gas ebullition was observed during ROV investigation [21, 63]. In the Shenhu area, acoustic blank in shallow horizons, presumably caused by the dissociation of deep-seated hydrates, is connected with the gas chimney [64, 65]. The derived AOM rate in Shenhu (17.5 mmol·m-2·yr-1) is consistent with a pioneer work (20.9 mmol·m-2·yr-1 and 11 mmol·m-2·yr-1) [27, 66]. In the Xisha area, gas chimneys, pockmarks, and gas bubble as well as chemosynthesis-based biota were also observed [32, 67]. These observations are consistent with the porewater profiles showing a bubble irrigation feature and shallow SMTZ in the Xisha area. The authigenic carbonate precipitation is due to the DIC production by AOM, resulting in the highest rate of carbonate precipitation in Xisha.Synthetic biogeochemical rates and fluxes of methane and DIC in the cold seeps and hydrate-bearing area of the northern SCS are also estimated by combining raster means of the three regions and arithmetical means of all the cores. In the northern SCS, area-based biogeochemical rates and fluxes of CH4 and DIC were far lower than estimates only based on few single points in the other regions (see Table 2). In the active continental margins, e.g., Hydrate Ridge, Costa Rica, and Hikurangi Margin, strong tectonic compression and structural fractures generally facilitate the ascension of internal methane-rich fluid across the sediment-seawater interface. Therefore, the rates of AOM and the fluxes of methane and DIC in the three areas are uniformly at least two orders of magnitude higher than those in the northern SCS [53, 58–60, 68, 69]. In the passive continental margins, e.g., Blake Ridge, Sakhalin Island, Skagerrak, and Ulleung Basin, AOM rates and fluxes of methane and DIC in the other regions worldwide are comparable with those in the northern SCS [37, 61, 70–72]. With the exception of the Gulf of Mexico and Vestnesa Ridge, the extraordinarily high AOM rates and the flux of methane can be attributed to cores retrieved near the gas chimney in the Gulf of Mexico [73] and around the pockmark in the Vestnesa Ridge [74]. As the world’s largest anoxic basin and the largest surface reservoir of aqueous methane, the Black Sea was known for its extremely high AOM rate and methane flux compared with the northern SCS [61]. Additionally, recent works have revealed that the AOM rate in marine sediments worldwide is ~7.9-10.7 mmol·m-2·yr-1 [75], which is less than the rate (39.9 mmol·m-2·yr-1) in the cold seep area of the northern South China Sea.Table 2
Areal comparison of biogeochemical rates and effluxes of methane and DIC with other cold seeps.
Location
CH4 efflux
DIC efflux
OSR
ME
AOM
CP
Reference
Northern SCS
5.7
46.2
16.5
3.4
39.9
13.7
This study
Hydrate Ridge
839.5
9000
51.83
6575
[60, 68, 69]
Gulf of Mexico
109.5
1718
[73]
Costa Rica
4562.5
5880
[59]
Black Sea
3017
11363
[61]
Blake Ridge
2.92
0.365
10.6
[37]
Sakhalin Island (Derugin Basin)
59.9
11.9
33.2
11
[37]
Skagerrak (Denmark)
36.5
62
[70]
Hikurangi Margin (New Zealand)
4380
5548
2.5
1.5
333
[53, 58]
Ulleung Basin
23
53.8
40
6.3
[71, 72]
Vestnesa Ridge
12.8
184.8
[74]


### 5.2. DIC and Aqueous Methane Turnover in the Shallow Sediments

Generally, the sources of DIC in cold seeps include OSR, AOM, and methanogenesis, whereas the major DIC sinks are carbonate precipitation and diffusive loss towards the bottom water [76, 77]. The relative contributions of the 5 processes mentioned above influencing the DIC pool in Dongsha, Xisha, Shenhu, and all cold seep sites in the northern SCS were shown in Figure 3. The three areas and the northern SCS exhibit similar patterns in terms of the sources and sinks of DIC pool. As for the sources of DIC, the contributions of each reaction to DIC production decrease in the order of AOM, OSR, and ME, with AOM accounting for more than 60% of the DIC production in all the three areas. This trend in the SCS is inconsistent with that in the Ulleung basin (AOM>ME>OSR), probably as a result of the difference in the reactivity of organic matter and modelled depth [71]. As for sinks of DIC, the fluxes of DIC towards the seafloor accounted for more than 80% of DIC sinks in all the areas and the authigenic carbonate precipitation only occupied about 20% of DIC sinks. Differences between sources and sinks in the three regions and northern SCS are 6.1 mmol·m-2·yr-1 in Dongsha, 17 mmol·m-2·yr-1 in Xisha, 2.7 mmol·m-2·yr-1 in Shenhu, and 5.9 mmol·m-2·yr-1 in the northern SCS. Even though mass conservations are kept in each modeling exercise, the imbalance in DIC may be attributed to the inclusion of numerous published calculations (see Table S5).Figure 3
Sources and sinks of methane and DIC in the shallow sediments. Positive values denote sources of methane or DIC and negative values denote sinks of methane or DIC. OSR: organoclastic sulfate reduction; AOM: anaerobic oxidation of methane; ME: methanogenesis; CP: carbonate precipitation; ES: external methane from the bubble dissolution within modeled sediment columns. ∆ represents difference between fluxes at the lower boundary and fluxes at the upper boundary. Percentages above or below each column represent the contributions of each biogeochemical reaction to the source/sink.E denotes the differences between sources and sinks of methane and DIC. Biogeochemical rates (OSR, AOM, ME, and CP) in the three regions and northern SCS are adopted from Table 1. The ES-CH4 was obtained by assuming a complete mass conservation of methane.There are four major factors, including methanogenesis, AOM, diffusive migration of methane, and phase transformation, that constrain methane cycling in shallow sediments. In Dongsha and Shenhu, AOM consumed a vast majority of methane moving up towards the seafloor. In contrast, AOM only consumed ~65% of methane and the residual was transported towards the seafloor in Xisha. As for the source of methane, bubble dissolution can widely take place if the dissolved methane concentration is lower than its phase equilibrium concentration [39]. In Xisha, bubble dissolution (99%) played a predominant role in methane supply and methanogenesis barely accounted for 1% of methane sources. However, constant diffusive supply of methane from the lower boundary provided the major source of methane in Dongsha and Shenhu. The diffusive migration of methane accounted for 47% and 82% of the methane source in Dongsha and Shenhu, respectively. Additionally, the methanogenesis played a similar role in supplying methane to the sediments (17% in Dongsha and 12% in Shenhu).


### 5.3. Influence of Shallow Carbon Cycling on the Carbon Reservoir of Seawater

As essential components of long-term carbon cycle, the confluence of organic carbon burial, organic matter degradation, and AOM mediates the carbon equilibrium between seawater reservoir and deep carbon and persistently affects the ocean environment [78]. Organic materials of both terrestrial and marine in origin undergo microbial decomposition during settling, and less than 1% of organic matter from the surface ocean reaches the seafloor [79]. Following the aerobic and anaerobic oxidation of organic matter in the sediments, the remaining refractory organic matter is permanently buried. According to the model results in each site, average POC burial flux below 20 mbsf is 186 mmol·m-2·yr-1 in the cold seeps and hydrate-bearing area of the northern SCS, which means that 2.8 Gmol POC are buried annually (see Figure 4). As the final step of the organic matter, methanogenesis produces both methane and some other dissolved organic matters (DOM), e.g., formate and acetate [80]. The majority of methane is consumed by AOM, thereby generating DIC as the byproduct. The remaining methane and other DOM may migrate towards the sediment-water interface [81].Figure 4
Shallow carbon cycling in the cold seeps and gas hydrate-bearing sediments in the northern South China Sea. The numbers in bold represent mean values of interpolation mean and arithmetic mean (Table1). The numbers in italic are derived from the model. The number in grey is estimated considering methane mass conservation.As we mentioned above, carbon entering the seawater from the sediments in the form of DIC and DOM can also, to some extent, influence the deep ocean carbon cycling, especially the carbonate saturation state in the water column [82, 83]. In the hydrate-bearing area of the northern SCS, DOM effluxes in the form of aqueous methane were ~5 mmol·m-2·yr-1 (see Figure 4). By multiplying by the area of 1.6×104km2, 9×10−5Tmol aqueous methane is released out of the sediments annually. Irrespective of gaseous methane in the intensive gas seepage sites, such considerable amounts of DOM diffusing towards the bottom water may play an important role in the DOM cycling in the deep waters of the northern SCS. Recent studies have shown much greater DOM fluxes (0.1 mol·m-2·yr-1) by accounting for other forms of DOM species such as volatile fatty acids [80, 83]. Thus, the widely distributed gas hydrates and cold seeps in the northern SCS could contribute far more than 9×10−5Tmol DOM/yr to the ocean if gaseous methane and other organic acids are incorporated into the flux calculation. Additionally, DIC efflux (46 mmol·m-2·yr-1) was several times higher than dissolved methane flux. Similarly, multiplying it by the area of cold seeps and hydrate-bearing areas gives an annual release amount of 7×10−4Tmol·C·yr−1. Despite that this estimation is not comparable to the DIC flux on the shelf [84], it may have an impact on the bottom water chemistry in cold seeps and hydrate-bearing areas. Release of DIC into bottom water can, in some case, prompt the production and preservation of biogenic and authigenic carbonate. Nevertheless, aerobic oxidation of methane released from sediments produces CO2 and lowers the seawater pH, thereby probably dissolving carbonate [82]. Accordingly, a total of 8.2×10−4Tmol dissolved carbon is released into bottom water annually but its influence on bottom water carbonate systems is still unclear.


## 5.1. Regional Biogeochemical Rates and Fluxes of Methane and DIC

By upscaling the rates and fluxes in the specific sites into a region of interest, raster means were obtained by ARCMAP based on the small subareas divided automatically, which represent areal means. To better constrain the areal mean values of individual biogeochemical reactions in each area, we averaged arithmetical mean and raster mean values (see Table1). Among the three areas, standard deviations of arithmetical means in the Xisha area were generally higher than 100, which indicate high heterogeneity of the biogeochemical rates. The arithmetical average of the rates of AOM in the Xisha area was much greater than those in the other two areas due to intensive methane bubbling. However, the raster mean of AOM rates in Xisha was slightly lower than those at the other two areas. The inconsistent result derived using two different methods is ascribed to the exclusion of intensive seepage sites due to unavailability of the sampling locations from the interpolation. Regardless of all the discrepancies between the two methods, a tentative comparison among the three regions indicates that the rates of AOM and carbonate precipitation in Xisha are overwhelmingly larger than those in the other two areas and that the rate of organic matter degradation in Xisha was slightly higher than those in Dongsha and Shenhu areas. The low rates of biogeochemical reactions in Shenhu are likely due to the absence of methane seepage and the occurrence of deep-seated hydrates [12]. The Dongsha area always shows the intermediate rates compared with the other two areas owing to the occurrence of sporadic methane seepage sites and gas hydrates in shallow sediments.Table 1
Biogeochemical rates and effluxes of methane and DIC in the three areas.
Location
CH4 efflux
DIC efflux
OSR
ME
AOM
CP
Cores
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
RM
AM
SD
Average
Dongsha
0.2
9.7
95
4.9
25.3
19.3
7
22.3
14.1
17.7
12
15.9
5.1
4.8
4
4.9
29.0
26.7
14
27.8
5.5
5.4
1
5.3
363
Xisha
0.1
144
366
72
36.0
246.2
136
141
27.6
14.6
20
21.1
2.6
1.5
1
2
9.4
261.1
385
135.2
27.3
37.4
28
32.3
8
Shenhu
0.0
0.1
1
0.0
17.1
19.8
5
18.4
5.8
9.8
2
7.8
1.4
2.8
1
2.1
16.8
18.2
5
17.5
5.4
6.3
1.6
5.8
24
North SCS
0.13
11.3
96
5.7
29.5
62.9
178
46.2
19.6
13.4
12
16.5
3.6
3.2
3
3.4
19.0
60.8
180
39.9
15.4
12
17
13.7
395
RM: raster mean; AM: arithmetical mean; SD: standard deviation.The Xisha area is characterized by high content of sedimentary organic matter (up to 2%) mainly sourced from marine algae [55, 56]. This may explain the relatively higher rate of organic matter degradation in Xisha compared to Dongsha and Shenhu. Applying the bimolecular kinetics to calculate the AOM rate, we assume that the adjustable rate constant (kAOM) represents the integrated effects of a number of factors including the microbial community structure and abundance, bioenergetics, and enzyme kinetics. Therefore, different settings may have distinct rate constants, ranging over 6 orders of magnitude [57]. Similar to the other passive continental margin environments, the rate constant (kAOM=1) in this study is generally lower than those in active continental margin environments [26, 37, 58–61]. In the cold seeps and hydrate-bearing area, methane expulsion is usually associated with some specific structures, e.g., mud diapir, fault, and gas chimney. These structures provide conduits for the upward migration of methane-rich fluid. According to the geophysical imaging in Dongsha, large amounts of mud diapir and gas chimney are well developed in subsurface sediments [62]. Gas bubble release was indicated by hydroacoustic flares, and indeed, gas ebullition was observed during ROV investigation [21, 63]. In the Shenhu area, acoustic blank in shallow horizons, presumably caused by the dissociation of deep-seated hydrates, is connected with the gas chimney [64, 65]. The derived AOM rate in Shenhu (17.5 mmol·m-2·yr-1) is consistent with a pioneer work (20.9 mmol·m-2·yr-1 and 11 mmol·m-2·yr-1) [27, 66]. In the Xisha area, gas chimneys, pockmarks, and gas bubble as well as chemosynthesis-based biota were also observed [32, 67]. These observations are consistent with the porewater profiles showing a bubble irrigation feature and shallow SMTZ in the Xisha area. The authigenic carbonate precipitation is due to the DIC production by AOM, resulting in the highest rate of carbonate precipitation in Xisha.Synthetic biogeochemical rates and fluxes of methane and DIC in the cold seeps and hydrate-bearing area of the northern SCS are also estimated by combining raster means of the three regions and arithmetical means of all the cores. In the northern SCS, area-based biogeochemical rates and fluxes of CH4 and DIC were far lower than estimates only based on few single points in the other regions (see Table 2). In the active continental margins, e.g., Hydrate Ridge, Costa Rica, and Hikurangi Margin, strong tectonic compression and structural fractures generally facilitate the ascension of internal methane-rich fluid across the sediment-seawater interface. Therefore, the rates of AOM and the fluxes of methane and DIC in the three areas are uniformly at least two orders of magnitude higher than those in the northern SCS [53, 58–60, 68, 69]. In the passive continental margins, e.g., Blake Ridge, Sakhalin Island, Skagerrak, and Ulleung Basin, AOM rates and fluxes of methane and DIC in the other regions worldwide are comparable with those in the northern SCS [37, 61, 70–72]. With the exception of the Gulf of Mexico and Vestnesa Ridge, the extraordinarily high AOM rates and the flux of methane can be attributed to cores retrieved near the gas chimney in the Gulf of Mexico [73] and around the pockmark in the Vestnesa Ridge [74]. As the world’s largest anoxic basin and the largest surface reservoir of aqueous methane, the Black Sea was known for its extremely high AOM rate and methane flux compared with the northern SCS [61]. Additionally, recent works have revealed that the AOM rate in marine sediments worldwide is ~7.9-10.7 mmol·m-2·yr-1 [75], which is less than the rate (39.9 mmol·m-2·yr-1) in the cold seep area of the northern South China Sea.Table 2
Areal comparison of biogeochemical rates and effluxes of methane and DIC with other cold seeps.
Location
CH4 efflux
DIC efflux
OSR
ME
AOM
CP
Reference
Northern SCS
5.7
46.2
16.5
3.4
39.9
13.7
This study
Hydrate Ridge
839.5
9000
51.83
6575
[60, 68, 69]
Gulf of Mexico
109.5
1718
[73]
Costa Rica
4562.5
5880
[59]
Black Sea
3017
11363
[61]
Blake Ridge
2.92
0.365
10.6
[37]
Sakhalin Island (Derugin Basin)
59.9
11.9
33.2
11
[37]
Skagerrak (Denmark)
36.5
62
[70]
Hikurangi Margin (New Zealand)
4380
5548
2.5
1.5
333
[53, 58]
Ulleung Basin
23
53.8
40
6.3
[71, 72]
Vestnesa Ridge
12.8
184.8
[74]


## 5.2. DIC and Aqueous Methane Turnover in the Shallow Sediments

Generally, the sources of DIC in cold seeps include OSR, AOM, and methanogenesis, whereas the major DIC sinks are carbonate precipitation and diffusive loss towards the bottom water [76, 77]. The relative contributions of the 5 processes mentioned above influencing the DIC pool in Dongsha, Xisha, Shenhu, and all cold seep sites in the northern SCS were shown in Figure 3. The three areas and the northern SCS exhibit similar patterns in terms of the sources and sinks of DIC pool. As for the sources of DIC, the contributions of each reaction to DIC production decrease in the order of AOM, OSR, and ME, with AOM accounting for more than 60% of the DIC production in all the three areas. This trend in the SCS is inconsistent with that in the Ulleung basin (AOM>ME>OSR), probably as a result of the difference in the reactivity of organic matter and modelled depth [71]. As for sinks of DIC, the fluxes of DIC towards the seafloor accounted for more than 80% of DIC sinks in all the areas and the authigenic carbonate precipitation only occupied about 20% of DIC sinks. Differences between sources and sinks in the three regions and northern SCS are 6.1 mmol·m-2·yr-1 in Dongsha, 17 mmol·m-2·yr-1 in Xisha, 2.7 mmol·m-2·yr-1 in Shenhu, and 5.9 mmol·m-2·yr-1 in the northern SCS. Even though mass conservations are kept in each modeling exercise, the imbalance in DIC may be attributed to the inclusion of numerous published calculations (see Table S5).Figure 3
Sources and sinks of methane and DIC in the shallow sediments. Positive values denote sources of methane or DIC and negative values denote sinks of methane or DIC. OSR: organoclastic sulfate reduction; AOM: anaerobic oxidation of methane; ME: methanogenesis; CP: carbonate precipitation; ES: external methane from the bubble dissolution within modeled sediment columns. ∆ represents difference between fluxes at the lower boundary and fluxes at the upper boundary. Percentages above or below each column represent the contributions of each biogeochemical reaction to the source/sink.E denotes the differences between sources and sinks of methane and DIC. Biogeochemical rates (OSR, AOM, ME, and CP) in the three regions and northern SCS are adopted from Table 1. The ES-CH4 was obtained by assuming a complete mass conservation of methane.There are four major factors, including methanogenesis, AOM, diffusive migration of methane, and phase transformation, that constrain methane cycling in shallow sediments. In Dongsha and Shenhu, AOM consumed a vast majority of methane moving up towards the seafloor. In contrast, AOM only consumed ~65% of methane and the residual was transported towards the seafloor in Xisha. As for the source of methane, bubble dissolution can widely take place if the dissolved methane concentration is lower than its phase equilibrium concentration [39]. In Xisha, bubble dissolution (99%) played a predominant role in methane supply and methanogenesis barely accounted for 1% of methane sources. However, constant diffusive supply of methane from the lower boundary provided the major source of methane in Dongsha and Shenhu. The diffusive migration of methane accounted for 47% and 82% of the methane source in Dongsha and Shenhu, respectively. Additionally, the methanogenesis played a similar role in supplying methane to the sediments (17% in Dongsha and 12% in Shenhu).


## 5.3. Influence of Shallow Carbon Cycling on the Carbon Reservoir of Seawater

As essential components of long-term carbon cycle, the confluence of organic carbon burial, organic matter degradation, and AOM mediates the carbon equilibrium between seawater reservoir and deep carbon and persistently affects the ocean environment [78]. Organic materials of both terrestrial and marine in origin undergo microbial decomposition during settling, and less than 1% of organic matter from the surface ocean reaches the seafloor [79]. Following the aerobic and anaerobic oxidation of organic matter in the sediments, the remaining refractory organic matter is permanently buried. According to the model results in each site, average POC burial flux below 20 mbsf is 186 mmol·m-2·yr-1 in the cold seeps and hydrate-bearing area of the northern SCS, which means that 2.8 Gmol POC are buried annually (see Figure 4). As the final step of the organic matter, methanogenesis produces both methane and some other dissolved organic matters (DOM), e.g., formate and acetate [80]. The majority of methane is consumed by AOM, thereby generating DIC as the byproduct. The remaining methane and other DOM may migrate towards the sediment-water interface [81].Figure 4
Shallow carbon cycling in the cold seeps and gas hydrate-bearing sediments in the northern South China Sea. The numbers in bold represent mean values of interpolation mean and arithmetic mean (Table1). The numbers in italic are derived from the model. The number in grey is estimated considering methane mass conservation.As we mentioned above, carbon entering the seawater from the sediments in the form of DIC and DOM can also, to some extent, influence the deep ocean carbon cycling, especially the carbonate saturation state in the water column [82, 83]. In the hydrate-bearing area of the northern SCS, DOM effluxes in the form of aqueous methane were ~5 mmol·m-2·yr-1 (see Figure 4). By multiplying by the area of 1.6×104km2, 9×10−5Tmol aqueous methane is released out of the sediments annually. Irrespective of gaseous methane in the intensive gas seepage sites, such considerable amounts of DOM diffusing towards the bottom water may play an important role in the DOM cycling in the deep waters of the northern SCS. Recent studies have shown much greater DOM fluxes (0.1 mol·m-2·yr-1) by accounting for other forms of DOM species such as volatile fatty acids [80, 83]. Thus, the widely distributed gas hydrates and cold seeps in the northern SCS could contribute far more than 9×10−5Tmol DOM/yr to the ocean if gaseous methane and other organic acids are incorporated into the flux calculation. Additionally, DIC efflux (46 mmol·m-2·yr-1) was several times higher than dissolved methane flux. Similarly, multiplying it by the area of cold seeps and hydrate-bearing areas gives an annual release amount of 7×10−4Tmol·C·yr−1. Despite that this estimation is not comparable to the DIC flux on the shelf [84], it may have an impact on the bottom water chemistry in cold seeps and hydrate-bearing areas. Release of DIC into bottom water can, in some case, prompt the production and preservation of biogenic and authigenic carbonate. Nevertheless, aerobic oxidation of methane released from sediments produces CO2 and lowers the seawater pH, thereby probably dissolving carbonate [82]. Accordingly, a total of 8.2×10−4Tmol dissolved carbon is released into bottom water annually but its influence on bottom water carbonate systems is still unclear.


## 6. Conclusion

This study is aimed at quantitatively assessing the carbon turnover within the shallow sediments and quantifying the areal effluxes of methane and DIC in Xisha, Dongsha, and Shenhu, northern SCS. A 1-D reaction-transport model was used to calculate the biogeochemical rates and fluxes of methane and DIC in the 37 new collected cores. In addition, a total of 395 cores, including 358 cores reported in pioneer literatures, were used to extrapolate the areal rates and fluxes by spatial interpolation in the northern SCS. The average of the arithmetical mean and interpolation mean revealed that the rates of AOM and carbonate precipitation and effluxes of methane and DIC in Xisha were at least one order of magnitude higher than those in Dongsha and Shenhu and the rate of organic matter degradation in Xisha was slightly higher than those in the other two areas. This is because of the occurrence of intensive methane gas bubbling in Xisha. However, the diffusive source played the predominant role in methane supply in the other two areas. The majority of DIC mainly produced by AOM is diffused out of sediments in the three areas. In comparison to those in the active continental margin and euxinic environment, the AOM rates and fluxes of methane and DIC in the cold seeps of the northern SCS were much lower. Alternatively, source and sink analysis of methane revealed that the bubble dissolution in Xisha contributed most of the methane. According to the estimates in this study, ~2.8 Gmol organic matter was buried and at least ~0.8 Gmol dissolved carbon was released from marine sediments annually. Nonetheless, to what extent the impact of carbon release from the shallow sediments on the oceanic carbon reservoir remains to be investigated.



---

*Source: 2573937-2019-03-17.xml*