A Quantitative Assessment of Methane-Derived Carbon Cycling at the Cold Seeps in the Northwestern South China Sea

: Widespread cold seeps along continental margins are significant sources of dissolved carbon to the ocean water. However, little is known about the methane turnovers and possible impact of seepage on the bottom seawater at the cold seeps in the South China Sea (SCS). We present seafloor observation and porewater data of six push cores, one piston core and three boreholes as well as fifteen bottom-water samples collected from four cold seep areas in the northwestern SCS. The depths of the sulfate–methane transition zone (SMTZ) are generally shallow, ranging from ~7 to < 0.5 mbsf (meters below seafloor). Reaction-transport modelling results show that methane dynamics were highly variable due to the transport and dissolution of ascending gas. Dissolved methane is predominantly consumed by anaerobic oxidation of methane (AOM) at the SMTZ and trapped by gas hydrate formation below it, with depth-integrated AOM rates ranging from 59.0 and 591 mmol m − 2 yr − 1 . The δ 13 C and ∆ 14 C values of bottom-water dissolved inorganic carbon (DIC) suggest discharge of 13 C- and 14 C-depleted fossil carbon to the bottom water at the cold seep areas. Based on a two-endmember estimate, cold seeps fluids likely contribute 16–26% of the bottom seawater DIC and may have an impact on the long-term deep-sea carbon cycle. Our results reveal the methane-related carbon inventories are highly heterogeneous in the cold seep systems, which are probably dependent on the distances of the sampling sites to the seepage center. To our knowledge, this is the first quantitative study on the contribution of cold seep fluids to the bottom-water carbon reservoir of the SCS, and might help to understand the dynamics and the environmental impact of hydrocarbon seep in the SCS. In present seafloor observations and porewater geochemical data of six push cores, one piston core and three boreholes as well as fifteen bottom-water samples collected at two newly discovered seeps with weak activities and an intensely active seep in the Qiongdongnan Basin. Using steady-state reaction-transport modeling, we aim at quantifying the methane turnover rates in shallow sediments at the seep sites. The δ 13 C DIC and Δ 14 C DIC values of bottom-water samples are also measured to investigate the potential contribution of fossil carbon released by cold seeps. Furthermore, a two-endmember simple mass balance model is applied to evaluate the contribution of cold seep fluids in this area. In combination with seafloor observations and geochemistry of fluids from cold seeps, we demonstrate the heterogeneity of fluxes methodology—J.F. and N.L.; formal analysis—J.F., N.L. and M.L.; resources—J.L., S.Y., H.W., D.C.; writing—original draft preparation —J.F. and N.L.; writing—review and editing —J.F. and M.L.; funding acquisition —S.Y., H.W., J.L. All authors


Introduction
The continental margin sediments contain large reservoirs of methane either as dissolved phase in porewater, solid gas hydrate or free gas (bubbles) depending on its in situ solubility. Submarine cold core and three boreholes as well as fifteen bottom-water samples collected at two newly discovered seeps with weak activities and an intensely active seep in the Qiongdongnan Basin. Using steady-state reaction-transport modeling, we aim at quantifying the methane turnover rates in shallow sediments at the seep sites. The δ 13 C DIC and ∆ 14 C DIC values of bottom-water samples are also measured to investigate the potential contribution of fossil carbon released by cold seeps. Furthermore, a two-endmember simple mass balance model is applied to evaluate the contribution of cold seep fluids in this area. In combination with seafloor observations and geochemistry of fluids from cold seeps, we demonstrate the heterogeneity of fluxes and geochemical processes at the cold seep areas of the Qiongdongnan Basin.

Geological Background
The northern SCS is characterized as a Cenozoic passive continental margin [32]. The Qiongdongnan Basin, located in the northwestern part of the SCS, is a northeastern trended Cenozoic sedimentary basin, which is covered by thick sedimentary materials [33]. This basin has experienced a rifting stage and a post-rift thermal subsidence stage. During the first stage, plenty of half-grabens and sags were formed in the basin. At the later stage, thermal subsidence occurred and thick sediment sequences mainly consisting of mudstones were deposited since Miocene. Consequently, the sedimentation rates and the geothermal gradients are both high in the Qiongdongnan Basin [34]. As a result, there are abundant hydrocarbon accumulations in this basin. Widespread faulting and/or diapirism facilitated the migration of hydrocarbon gas to the upper strata [35]. Gas hydrate reservoirs were identified via numerous bottom-simulating reflectors (BSRs) and gas chimneys as the main types of fluid conduit in this region [36][37][38].
The Haima cold seeps have been discovered in the southern uplift belt of the Qiongdongnan Basin with water depths of~1400 m during R/V Haiyang-6 cruises in 2015 and 2016. Abundant chemosynthetic communities, methane-derived authigenic carbonates and massive gas hydrates in the near-surface sediments were observed within the seepage area [23]. Seafloor observations and sampling were conducted at sites HM-ROV05 and HM-ROV within the Haima cold seeps. In addition, another two weak seeps (QH-ROV05 and QH-ROV07) were also dived during R/V Haiyang-6 cruises in 2018 ( Figure 1). These seeps are located above two acoustic blanking zones that were inferred as gas chimneys to the northeast of the Haima cold seeps (Figure 1). The bathymetry in these areas is characterized by flat topography with water depths ranging from 1700 to 1800 m (Figure 1). During a later gas hydrate drilling expedition, multiple gas hydrates were observed from 8 to 174 mbsf of site QH-W08-2018 (hereinafter referred to as site W08, including two nearby holes W08B and W08C) within the investigation area QH-ROV07 and from 15 to 160 mbsf of site QH-W09-2018 (hereinafter referred to as W09) within the area QH-ROV05, respectively [28,39,40].
Minerals 2020, 10, x FOR PEER REVIEW 3 of 25 and the discharge of considerable amount of methane and DIC to the water column around the cold seeps [29][30][31]. In this study, we present seafloor observations and porewater geochemical data of six push cores, one piston core and three boreholes as well as fifteen bottom-water samples collected at two newly discovered seeps with weak activities and an intensely active seep in the Qiongdongnan Basin. Using steady-state reaction-transport modeling, we aim at quantifying the methane turnover rates in shallow sediments at the seep sites. The δ 13 CDIC and Δ 14 CDIC values of bottom-water samples are also measured to investigate the potential contribution of fossil carbon released by cold seeps. Furthermore, a two-endmember simple mass balance model is applied to evaluate the contribution of cold seep fluids in this area. In combination with seafloor observations and geochemistry of fluids from cold seeps, we demonstrate the heterogeneity of fluxes and geochemical processes at the cold seep areas of the Qiongdongnan Basin.

Geological Background
The northern SCS is characterized as a Cenozoic passive continental margin [32]. The Qiongdongnan Basin, located in the northwestern part of the SCS, is a northeastern trended Cenozoic sedimentary basin, which is covered by thick sedimentary materials [33]. This basin has experienced a rifting stage and a post-rift thermal subsidence stage. During the first stage, plenty of half-grabens and sags were formed in the basin. At the later stage, thermal subsidence occurred and thick sediment sequences mainly consisting of mudstones were deposited since Miocene. Consequently, the sedimentation rates and the geothermal gradients are both high in the Qiongdongnan Basin [34]. As a result, there are abundant hydrocarbon accumulations in this basin. Widespread faulting and/or diapirism facilitated the migration of hydrocarbon gas to the upper strata [35]. Gas hydrate reservoirs were identified via numerous bottom-simulating reflectors (BSRs) and gas chimneys as the main types of fluid conduit in this region [36][37][38].
The Haima cold seeps have been discovered in the southern uplift belt of the Qiongdongnan Basin with water depths of ~1400 m during R/V Haiyang-6 cruises in 2015 and 2016. Abundant chemosynthetic communities, methane-derived authigenic carbonates and massive gas hydrates in the near-surface sediments were observed within the seepage area [23]. Seafloor observations and sampling were conducted at sites HM-ROV05 and HM-ROV within the Haima cold seeps. In addition, another two weak seeps (QH-ROV05 and QH-ROV07) were also dived during R/V Haiyang-6 cruises in 2018 ( Figure 1). These seeps are located above two acoustic blanking zones that were inferred as gas chimneys to the northeast of the Haima cold seeps ( Figure 1). The bathymetry in these areas is characterized by flat topography with water depths ranging from 1700 to 1800 m ( Figure 1). During a later gas hydrate drilling expedition, multiple gas hydrates were observed from 8 to 174 mbsf of site QH-W08-2018 (hereinafter referred to as site W08, including two nearby holes W08B and W08C) within the investigation area QH-ROV07 and from 15 to 160 mbsf of site QH-W09-2018 (hereinafter referred to as W09) within the area QH-ROV05, respectively [28,39,40].

Seafloor Observations
Two remotely-operated vehicle (ROV) dives were conducted around sites W08 and W09 by ROV "Haima" in April 2018 and Fugro ROV "FCV 2000D" in June 2018. Another dive was conducted in May 2018 at the seep site C using the ROV "Haima". Photos of authigenic carbonates and chemosynthetic communities and/or gas bubbles at these sites were taken during these cruises. Moreover, methane concentrations in bottom seawater were measured using a methane sensor mounded on the ROV "Haima".

Sampling and Analytical Methods
During the R/V Haiyang-6 cruise, 6 push cores (~70 cm) and 15 bottom-water samples were collected from sites QH-ROV05, QH-ROV07, HM-ROV05 and HM-ROV using the ROV "Haima" in 2018 ( Figure 1, Table 1). In QH-ROV07, 2 push cores, R7 and R7-2, were recovered at the edge of an authigenic carbonate deposit, where the gas hydrate-bearing core W08C was drilled only~5 m away. Another hydrate-bearing core W08B was drilled~20 m away from the carbonate deposit where a dome-like structure and fragments of recently dead seep-associated bivalves were observed. Two other push cores, R7-1 and R7-3, were taken~40 m away from W08B. In QH-ROV05, a push core ROV05-1 and a piston core (CL48) were also collected about 200 m and 1.1 km away from site W09, respectively ( Figure 1). In HM-ROV, a push core HM-01 was collected near the seep C with acoustic flare (Figure 1). No sediment core was collected at site HM-ROV05 due to a lack of time during that cruise. On the other hand, bottom-water samples, including R01-2018, ROV07, R-07, R-07-1 around W08, R-05-shell and ROV05 around W09, HM-R003-1 from HM-ROV05, HM-2-vent and HM-3-vent from HM-ROV were collected using specially-made fluid samplers which are temperature-held and pressure-tight~5 m above the seafloor. Other bottom-water samples, including ROV05-1 from W09, HM-1 from HM-ROV05 as well as HM-2, HM-3 from HM-ROV, were collected by Niskin bottles~5 m above the seafloor. In addition, three other samples, including ROV07+v, ROV07-1 from W08 and ROV05-1 from W09, were collected from the top of push cores ( Table 1).
The recovered cores were immediately brought to the onboard laboratory for porewater extraction. The porewater samples of the push cores and the piston core were collected at 5 or 10 cm intervals and at 60 cm intervals using Rhizon samplers, respectively. Porewater samples were acidified with ultra-pure concentrated HNO 3 (10 µL HNO 3 per 1 mL sample) for determining the concentrations of major elements onshore. Porewater samples for analysis of the concentration and carbon isotopic composition of DIC were preserved with a saturated HgCl 2 solution (∼20 µL HgCl 2 per 5 mL of sample). All the porewater samples were stored at 4 • C until further analysis. In addition, the porewater PO 4 3− concentrations of the piston core CL48 were measured onboard using the spectrophotometric method according to [41] with a UV-Vis Spectrophotometer (Hitachi U5100, Hitachi Limited, Tokyo, Japan). The precision for phosphate was ±3.0%. The total alkalinity (TA) of CL48 was determined onboard by direct titration with~0.006 M HCl using a pH meter. The analysis was calibrated using the seawater standard of International Association for the Physical Sciences of the Oceans (IAPSO), with a precision and detection limit of 0.05 meq L −1 . For the bottom-water samples, total alkalinity (TA) contents were determined onboard by the same method.
Moreover, porewater sampling and analysis of sulfate (SO 4 2− ) concentrations in the porewater samples of W08 and W09 were reported in [28]. Calcium (Ca 2+ ) concentrations were also measured aboard via ion chromatography (Dionex ICS-2100, Thermo Fisher Scientific, Waltham, MA, USA). The 4 ml sediment samples of W08 and W09 were sealed in 26 mL glass vials onboard. The vials were placed in a 60 • C oven for 30 min to fully release the hydrocarbon gases weakly adsorbed on the surfaces and pores of the sediment particles. After that, the concentrations of hydrocarbon gas in the headspace vials were measured onboard using the gas chromatograph method (Inficon Fusion MicroGC). The detection limit for all gases was 5 ppm. For the porewater samples from the push and piston cores, SO 4 2− and Ca 2+ concentrations were measured on a Dionex ICS-5000+ ion chromatograph (analytical precision of <2%) at the South China Sea Institute of Oceanology, Chinese Academy of Sciences. The anion (SO 4 2− ) and cation (Ca 2+ ) concentrations were determined by 500-and 100-fold dilution, respectively, using ultra-pure water. Concentration and isotopic analyses of DIC were determined using a Thermo Finnigan Gas Bench coupled with a Thermo Finnigan Delta V Advantage at the Louisiana State University. The analytical precisions were better than 0.1% for δ 13 C. For the porewater samples from W08 and W09, DIC concentrations and δ 13 C DIC values were determined via a continuous flow mass spectrometer (Thermo Delta V Advantage). The analytical precisions were better than 0.2% for δ 13 C. All δ 13 C data in this study are reported in per mil (% ) using the δ notation, relative to the standards Vienna-Pee Dee Belemnite (V-PDB). The natural radiocarbon 14 C contents of DIC (∆ 14 C DIC ) for bottom-water samples were obtained at Beta Analytic Inc., Miami, FL, USA, using accelerator mass spectrometry (AMS). Sample pretreatment, preparation, and measurement were conducted at Beta Analytic Inc. The Cambridge half-life (5730 ± 40 years) was used to calculate the apparent radiocarbon age and ∆ 14 C [42,43]. The accuracy and precision of the ∆ 14 C DIC measurements was 0.1% and 1.7% , respectively.

Reaction-Transport Model
In this study, a one-dimensional, steady-state reaction-transport model was applied to simulate one solid (POC) and four dissolved species including sulfate (SO 4 2− ), methane (CH 4 ), dissolved inorganic carbon (DIC) and calcium (Ca 2+ ). The model is modified from previous simulations of methane-rich sediments [20,[44][45][46], and a full description of the model is shown in the Appendix A.
All the reactions considered in the model and the expressions of kinetic rate are listed in Table 2.
is the POC degradation rate, a 0 (yr) is the initial age of organic matter in surface sediments, ν s (cm yr −1 ) is the burial velocity of solids, x (cm) is the depth in the sediment, Kc is an inhibition constant for POC degradation due to DIC and CH 4 accumulation in the porewater, POC (wt.%) is the POC content in sediments. R SR (mmol cm −3 yr −1 of SO 4 2− ) is the rate of sulfate reduction, [SO 4 2− ] is the SO 4 2− concentration, K SO 2 4 is the Michaelis-Menten constant for the inhibition of sulfate reduction at low sulfate concentrations, f POC converts between POC (dry wt.%) and DIC (mmol cm −3 of porewater): f POC = MW C /10Φ/(1 − Φ)/ρ S , where MWC is the molecular weight of carbon (12 g mol −1 ), ρ S is the density of dry sediments, and Φ is the porosity. R MG (mmol cm −3 yr −1 of CH 4 ) is the rate of methanogenesis. R AOM (mmol cm −3 yr −1 of CH 4 ) is the rate of AOM, k AOM is the rate constant, [CH 4 ] is the concentration of dissolved CH 4 . R CP (mmol cm −3 yr −1 of Ca 2+ ) is the rate of authigenic carbonate precipitation, k Ca (mol·cm −3 ·yr −1 ) is the rate constant, K SP (mol 2 ·L 2 ) is the thermodynamic equilibrium constant, [Ca 2+ ] and [CO 3 2− ] are the concentrations of Ca 2+ and CO 3 2− , respectively.
Due to the porewater sampling resolution, any influence of bioturbation and bioirrigation on the porewater profiles cannot be resolved. Therefore, solutes are assumed to be mainly transported by molecular diffusion and porewater burial, while solid species are assumed to be transported only by burial with prescribed compaction. At sites W08C and W09, the porewater species concentrations were close to that of seawater in the upper~2 and~3 m, respectively. This feature may be attributed to the bubble irrigation which is often observed at cold seep sites [20,30,45]. The porewater mixing with bottom water induced by rising gas bubbles can be described as a nonlocal transport similar to bioirrigation.
where α 0 (yr −1 ) is the coefficient of irrigation intensity, L irr (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 L irr , C 0 is the solute concentration at the SWI, and C x is the concentration at any depth within the irrigation zone.
Major biogeochemical reactions considered in the model are particulate organic matter (POM) degradation via sulfate reduction, AOM, methanogenesis and authigenic carbonate precipitation. Organic matter mineralization via aerobic respiration, denitrification, and metal oxide reduction were ignored because these processes mainly occur in the surface centimeter-scale sediments, which cannot be resolved by our sampling resolution.
The total rate of POM mineralization, R POC (wt.% C yr −1 ), is calculated via the power law model in which the initial age of organic matter in surface sediments, a 0 (yr) is considered in surface sediments [47]. a 0 is constrained using the measured PO 4 3− concentrations of the reference core.
When sulfate is almost completely consumed, the remaining POM is degraded to CO 2 and CH 4 via methanogenesis: The main pathways of methanogenesis in marine sediments are organic matter fermentation and CO 2 reduction [48]. Their net reactions at steady state are balanced as equivalent amounts of CO 2 and CH 4 are produced when per mole of POM is degraded [49]. Accordingly, the reaction of methanogenesis is a net reaction.
Methane is considered to be consumed by AOM [3]: The rate constant for AOM, k AOM (cm 3 mmol −1 yr −1 ), is tuned to the sulfate profiles within the SMTZ.
The loss of Ca 2+ resulting from the precipitation of authigenic carbonates as calcite (Ca 2+ + HCO 3 − → CaCO 3 + H + ) was simulated in the model using the thermodynamic solubility constant as defined by [50] (Table 2). A porewater pH value of 7.3 was used to calculate CO 3 2− from modeled DIC concentrations [51]. CaCO 3 was not simulated explicitly in the model. The length of the simulated model domain was set to 1000 cm for W08B and W08C, 2000 cm for W09, 500 cm for R7-1 and R7-3 and HM-1. Upper boundary conditions for all species were imposed as fixed concentrations (Dirichlet boundary) using measured values in the uppermost sediment layer where available. A zero concentration gradient was imposed at the lower boundary for all the species except CH 4 . CH 4 concentration at the lower boundary was a tunable parameter constrained from the SO 4 2− profile. The model was solved using the NDSolve object of MATHEMATICA V. 10.0 (Wolfram Research, Champaign, IL, USA). The steady-state simulations were run for 10 7 years to achieve the steady state with a mass conservation of >99 %. Further details on the model solutions can be found in Appendix A.

Site Characteristics
Seafloor observations showed that a small mud mound sparsely colonized by living clams at site QH-ROV05 (Figure 2). At site QH-ROV07, massive authigenic carbonate deposits and dead clams were observed on the seabed ( Figure 2B,C). Gas hydrate-bearing core W08B was drilled~20 m away from the carbonate deposit where there is a dome-like structure ( Figure 2D). Fragments of recently dead seep-associated bivalves and tubeworms ( Figure 2D) were scattered on the periphery of the dome. Scatter seep-associated bivalve fragments were observed on the seafloor. The sediments of the cores R7-1 and R7-3 mainly consist of black-green silty clay with some shell fragments and a strong odor of hydrogen sulfide.

General Geochemical Trends
Profiles of porewater SO4 2− , DIC, Ca 2+ concentrations are shown in Figures 3-5, Tables B1 and B2. In the borehole W08B, the SO4 2− concentrations were nearly depleted near the seafloor, whereas DIC concentrations were high and Ca 2+ concentrations were depleted from the seabed below. In the borehole W08C, the SO4 2− concentrations were depleted below 8 mbsf, which indicated the depth of SMTZ was shallow and located at ~2-8 mbsf. Due to the coarse sampling resolution in the upper 10 At site HM-ROV05, there were extensive authigenic carbonate pavements with bacterial mats, squat lobsters and a small amount of living tubeworms colonizing the fractures ( Figure 2E). Dead bivalves were scattered and gas bubbling was observed on the mud mounds ( Figure 2F,G). Abundant deep-sea mussels were clustered on the flank of the mud mounds ( Figure 2H). Methane concentrations in bottom water were <5.4 nM at site QH-ROV07, 5.4-36 nM at site QH-ROV05, 5.4-900 nM at site HM-ROV05 and >5.4 × 10 4 nM at site HM-ROV, respectively.
In the borehole W08B, the SO 4 2− concentrations were nearly depleted near the seafloor, whereas DIC concentrations were high and Ca 2+ concentrations were depleted from the seabed below. In the borehole W08C, the SO 4 2− concentrations were depleted below 8 mbsf, which indicated the depth of SMTZ was shallow and located at~2-8 mbsf. Due to the coarse sampling resolution in the upper 10 m of the hydrate-bearing borehole W08C, no clear downcore trend of porewater geochemical data was observed (Table A3). Considering the close proximity between the push cores (R7 and R7-2) and W08C, we present the porewater data of these cores in one panel for the sake of modelling. In contrast, both R7 and R7-2 do not show significant downcore variations in SO 4 2− , DIC and Ca 2+ concentrations.   Table A3). The profiles δ 13 C DIC values mirrored those of DIC concentrations reaching minimum at the SMTZ of the cores and shifting to positive values in the methanogenic zone as shown in the porewater profiles of W08B, W08C and W09 (Figures 3-5, Table A3). The δ 13 C DIC values of R7 and R7-2 showed slightly downcore decrease, whereas those of R05-1 and R05-2 display near-seawater values with depth ( Figures 3-5, Table A4). Natural radiocarbon measurements of DIC for bottom-water samples yielded the 14 C ages of DIC ranging from 1250 to 590 years BP (∆ 14 C = −151% to −71% ). The ∆ 14 C DIC values of temperatureand pressure-tight (T,P-tight) samples range from −149% to −95% (Table A5). Moreover, the δ 13 C DIC values of T,P-tight samples (−3.4% to −1.6% ) were generally lower than those of T,P-tight free samples (−2.1% to −1.7% ). Small variations were measured in pH values and total alkalinity (TA) concentrations of the bottom-water samples (values ranged from 7.6 to 7.9 for pH and ranged from 2.8 to 3.2 mM for TA) (Table A5).

Reaction-Transport Modeling
The simulation profiles and methane turnover rates are shown in Figures 3-5 and Table 3, respectively. The model parameters are listed in Table A2. The porewater geochemistry of the hydrate-bearing hole W08B displayed that the SMTZ must be located near the seabed with depth shallower than 1 mbsf, but its exact depth cannot be resolved by the sampling scheme during the gas hydrate drilling expedition (Figure 3, Table 3). The CH 4 and SO 4 2− profiles and methane turnover were instead constrained by the depth where gas hydrates first occurred and may represent the minimum values. Our simulations generally reproduced the measured concentrations of SO 4 2− , DIC and Ca 2+ in the investigated cores except the enigmatic reversals in R7-3 below 0.5 mbsf.
Minerals 2020, 10, x FOR PEER REVIEW 11 of 25  The rates of POC degradation through sulfate reduction ranged between 0.1 and 1.8 mmol m −2 yr −1 at the study cores. Compared to the low depth-integrated rates of POC degradation, the AOM overwhelmingly dominated the depletion of sulfate with depth-integrated rates of 591, 383, 241, 410, 80.3 and 59.0 mmol m −2 yr −1 for cores W08B, W08C, W09, HM-1, R7-1 and R7-3, respectively. The AOM rates were mainly sustained by an external methane source, and methanogenesis contributed only a negligible amount of methane (Table 3). Noting that the methane turnovers of W08B may be underestimated due to the lack of data in the surface sediments and the uncertainty of the exact depth of SMTZ. The benthic DIC fluxes were estimated to be 460, 272, 211, 295, 57.1 and 50.8 mmol m −2 yr −1 for cores W08B, W08C, W09, HM-1, R7-1 and R7-3, respectively (Table 3). Natural radiocarbon measurements of DIC for bottom-water samples yielded the 14 C ages of DIC ranging from 1250 to 590 years BP (Δ 14 C = −151‰ to −71‰). The Δ 14 CDIC values of temperature-and pressure-tight (T,P-tight) samples range from −149‰ to −95‰ (Table B3). Moreover, the δ 13 CDIC values of T,P-tight samples (−3.4‰ to −1.6‰) were generally lower than those of T,P-tight free samples (−2.1‰ to −1.7‰). Small variations were measured in pH values and total alkalinity (TA) concentrations of the bottom-water samples (values ranged from 7.6 to 7.9 for pH and ranged from 2.8 to 3.2 mM for TA) (Table B3).

Reaction-Transport Modeling
The simulation profiles and methane turnover rates are shown in Figures 3-5 and Table 3, respectively. The model parameters are listed in Table A2. The porewater geochemistry of the hydrate-bearing hole W08B displayed that the SMTZ must be located near the seabed with depth shallower than 1 mbsf, but its exact depth cannot be resolved by the sampling scheme during the gas hydrate drilling expedition ( Figure 3, Table 3). The CH4 and SO4 2− profiles and methane turnover were instead constrained by the depth where gas hydrates first occurred and may represent the minimum values. Our simulations generally reproduced the measured concentrations of SO4 2− , DIC and Ca 2+ in the investigated cores except the enigmatic reversals in R7-3 below 0.5 mbsf.
The rates of POC degradation through sulfate reduction ranged between 0.1 and 1.8 mmol m −2 yr −1 at the study cores. Compared to the low depth-integrated rates of POC degradation, the AOM overwhelmingly dominated the depletion of sulfate with depth-integrated rates of 591, 383, 241, 410, 80.3 and 59.0 mmol m −2 yr −1 for cores W08B, W08C, W09, HM-1, R7-1 and R7-3, respectively. The AOM rates were mainly sustained by an external methane source, and methanogenesis contributed only a negligible amount of methane ( Table 3). Noting that the methane turnovers of W08B may be underestimated due to the lack of data in the surface sediments and the uncertainty of the exact depth of SMTZ. The benthic DIC fluxes were estimated to be 460, 272, 211, 295, 57.1 and 50.8 mmol m -2 yr -1 for cores W08B, W08C, W09, HM-1, R7-1 and R7-3, respectively (Table 3).

Methane-Related Carbon Cycling at Cold Seep Areas
In the Qiongdongnan Basin, upward expulsion of free gas in the sediment was identified by the blanking or pull-up seismic reflections in gas chimney and mud diapir structures [26,37,38]. Previous studies have shown that hydrocarbon seeps of W08 and W09 are characterized by abundant thermogenic gas within gas chimneys and a lack of advective fluid flow [28,39,40]. On the other hand, biogenic gas was transported upwards within fluid conduits and emitted to the water column at the stations HM-ROV and ROV2 in the eastern part of the Haima cold seeps [23,26]. The biogenic origin of emitted methane is also evident by extremely-depleted 13 C DIC in porewater profiles of core HM-01 with the lowest δ 13 C DIC of~−50% ( Figure 5; Table A4). In addition, the δ 13 C DIC values below the SMTZ became more positive (Figures 3-5), which is caused by the generation of 13 C-enriched DIC via methanogenesis [52,53]. Nevertheless, modeling results display that in-situ methanogenesis rates in the upper sediments are too low to supply sufficient methane to form gas hydrate (Table 3). Hence, thermogenic or microbial gas transport from deep-seated sediments serves as the main methane source for AOM and hydrate formation in the shallow sediments of the investigated cores.
Our modeling results also show that, compared to sulfate reduction via POC degradation (OSR) throughout the sulfate reduction zone, AOM at the bottom of the sulfate reduction zone is the major pathway for the consumption of dissolved sulfate ( Table 4). Noting that the methane-related carbon inventories are highly heterogeneous within short distances in an individual cold seep system. Outside of the seepage center within the gas chimney of QH-ROV07 (cores R7-1, R7-3), QH-ROV05 (CL48), HM-ROV (QDN-31) and HM-ROV05 (QDN-14A, QDN-14B), almost all the dissolved gas is consumed by AOM at lower rates, or there are no apparent upward flux of dissolved gas (Table 3; Table 4). In contrast, for the cores (W08B and W08C; W09; HM-01) within the seepage center, AOM rates are at least one order-of magnitude higher (Table 3; Table 4). Modeling results indicate that, above the GHOZ, AOM within the SMTZ is the main methane sink and also consume the majority of methane at these cores (Table 3). AOM increased porewater alkalinity by producing bicarbonate and led to authigenic carbonate precipitation as indicated by the decrease in Ca 2+ concentrations with depth (Figures 4-6). Notation: F SO4 and F CH4 are the downward flux of SO 4 2− and the upward flux of CH 4 , while R AOM refers to the depth-integrated reaction rate of AOM (unit: mmol m −2 yr −1 ). Z SMTZ is the depth of SMTZ.
Comparing our model results with those of other sites in the Qiongdongnan Basin and in other areas of the SCS, the methane turnover rates in shallow sediments are much higher at the cold seeps in the Qiongdongnan Basin than those at a dormant pockmark in the nearby SW Xisha Uplift, as well as those at sites located on deep-seated gas hydrate reservoirs in the Shenhu area and the Beikang Basin (Table 4). Besides, the methane turnover rates at the active cold seep sites in the Taixinan Basin are of the same order of magnitude as those in the Qiongdongnan Basin (Table 4). The porewater profiles in the shallow sediments in these two basins often showed kink-type shape, which are attributed to irrigation of gas bubbling or recent increase in upward methane flux related to non-steady states of fluid seepage [20,30,45]. Considering the geological settings of the sampling sites, it is suggested that differences in the proximity of the sites to the fluid conduits including fractures, faults or anticlines mainly account for the variability in upward methane fluxes and turnover rates [20,28,31,39,45]. Overall, combing seafloor observations with geochemical modeling, our results further demonstrate that the methane-related carbon inventories are highly heterogeneous at different cold seep systems or different parts within short distances in an individual cold seep system. This signature should be attributed to the difference in the distance between the sampling sites and the seepage center.

Potential Contribution of Fossil Carbon from Cold Seeps to Bottom-Water Carbon Pool
Cold seep systems contribute considerable proportions of the local carbon budget of the overlying water column by emitting large amount of fossil carbon with depleted 13 C and 14 C into the water column. The seepage often results in a significant decrease in the δ 13 C DIC and ∆ 14 C DIC , and a small increase in the DIC in the overlying seawater, either by in-situ oxidation of vent methane or the concurrent input of DIC from seeps, or both [10,[58][59][60][61]. Generally, the study of the DIC system in seawater can be greatly simplified by describing the system in terms of the alkalinity. In shallow sediments where pH is between 7.1 and 8.1, total alkalinity is often treated as carbonate alkalinity by ignoring other minor species for practical purposes [62].  [63], whereas the shade in B represents the background Δ 14 CDIC values of seawater of SCS [64]. The small increase in DIC in the seep sites of the Gulf of Mexico (GOM) and Cascadia convergent margin was attributed to the removal of DIC by active deposition of authigenic carbonates [10,60]. The dataset includes those of the bottom water (and overlying water column) and/or porewater at the seep sites of the Okinawa Trough [65], GOM [10,66], Hudson Canyon [61], Cascadia convergent margin [60], Western Sagami Bay [58,59], Santa Monica Basin [67] and Bullseye vent [68], as well as the seep-free sites of the NSCS [63].
To understand the impact of a seep-derived source of TA (excess alkalinity) on deep water inorganic carbon pool, the following simple two-endmember mass balance model is used: where the subscript "wc", "ex" and "bg" stand for the water column, excess TA, and background value of the SCS, respectively. The background concentration and δ 13 C of TA in the bottom water of  [63], whereas the shade in B represents the background ∆ 14 C DIC values of seawater of SCS [64]. The small increase in DIC in the seep sites of the Gulf of Mexico (GOM) and Cascadia convergent margin was attributed to the removal of DIC by active deposition of authigenic carbonates [10,60]. The dataset includes those of the bottom water (and overlying water column) and/or porewater at the seep sites of the Okinawa Trough [65], GOM [10,66], Hudson Canyon [61], Cascadia convergent margin [60], Western Sagami Bay [58,59], Santa Monica Basin [67] and Bullseye vent [68], as well as the seep-free sites of the NSCS [63].
Studies show that the benthic DIC fluxes at the SWI are in the order of magnitude of 10-10 3 mmol m −2 yr −1 in the cold seep areas of the Qiongdongnan Basin (Table 4). In addition, compared with other seep sites, the slightly elevated alkalinity concentrations together with negative δ 13 C values of bottom water suggest cold seep system may have influenced the carbon pool of bottom water ( Figure 6, Table 4). In addition, the depletion in 14 C of DIC (−178% to −71% ) suggests that oxidation of CH 4 from deeper reservoirs is likely the source of ancient carbon to bottom waters ( Figure 6B, Table A5).
To understand the impact of a seep-derived source of TA (excess alkalinity) on deep water inorganic carbon pool, the following simple two-endmember mass balance model is used: TA wc δ 13 C wc = TA ex δ 13 C ex + TA bg δ 13 C bg , where the subscript "wc", "ex" and "bg" stand for the water column, excess TA, and background value of the SCS, respectively. The background concentration and δ 13 C of TA in the bottom water of northern SCS are ca. 2.4 mM and 0 % , respectively [63]. By inserting the background and our measured values of TA into the above two equations, we estimate the TA ex and its carbon isotope composition are approximately 0.5-0.8 mM and −22.7% to −4.0% , respectively ( Figure 6; Table 4). The δ 13 C compositions of the excess carbonate alkalinity are similar to those of the pore fluids at the top of the push core near the cold seep vent (−19.8% at HM-1, Table A3). The excess δ 13 C values of DIC are to be comparable to those of vent fluids or uppermost pore fluids on the Cascadia convergent margin (−15.6% to −3.9% ) [60], on the slope of Gulf of Mexico (−27.8% ) [10], and in the Guaymas Basin (−25.6% ) [69]. These characteristics support the hypothesis that the carbonate alkalinity and δ 13 C DIC anomalies in the water samples result from the mixture of seawater and cold seep fluids emitting on the seafloor. Moreover, the bottom-water samples with the lowest δ l3 C DIC values were collected in the stations with living or recently dead chemosynthetic bivalves, further indicating that high-flux seep sites serve as an important source of DIC to the water column ( Figure 6A, Table A5).Model results show that benthic DIC fluxes in the study area are 10 1 -10 3 mmol m −2 a −1 , which is in the range of those reported at other cold seep sites in the SCS [22,30,45,46,56,57]. We postulate that the DIC released from seep-impacted sediments could alter the dissolved carbon pool in the overlying bottom water. To further quantify the local DIC contributions from cold seeps in the study area, we estimate the ratio of the excess alkalinity to the alkalinity of the bottom-water samples. The calculations yield the contribution of DIC from cold seep fluids to the total bottom water ranging between 16% and 26% (Table A5). This is similar to the estimates in the Gulf of Mexico and Okinawa Trough where cold seep fluids can contribute up to 14.3% and 20% of excess DIC to the bottom water, respectively [10,65]. Hydroacoustic imaging revealed that the height of the gas plume was~770 m at seep C within the Haima cold seeps [26], yet the contribution of cold seep fluids to the water column is still unknown. Here our calculation suggests that the upward DIC flux from the cold seeps could contribute up to one quarter of 13 C-depleted DIC to the bottom-water inorganic carbon pool. Therefore, the impact of cold seep fluids to the overlying water column inorganic carbon reservoir seems to be significant.
Previous investigation has showed that the area of seepage at the Haima cold seeps is approximately 350 km 2 [70]. Based on the cores sampled from this area by far, the model-derived benthic DIC fluxes in this area ranged between 15 and 1,139 mmol m −2 yr −1 with an average of 239 mmol m −2 yr −1 (Table 4). By multiplying this value by the area of seepage, the average DIC efflux amounts to~8.4 × 10 −5 Tmol yr −1 at this active cold seeps. The amount of DIC derived from the cold seeps is likely greater than our estimation because the DIC efflux is highly heterogeneous in a cold seep area. It is suggested that release of DIC into bottom water can, in some cases, promote the production and preservation of biogenic and authigenic carbonate. Nevertheless, aerobic oxidation of methane from subseafloor can produce CO 2 and lower the seawater pH, thereby probably dissolving carbonate [71]. A previous study showed that, at least~7 × 10 −4 Tmol DIC was released from marine sediments per year in cold seeps and hydrate-bearing areas assuming an area of 1.6 × 10 4 km 2 in the northern SCS [22]. Compared to this estimation, our results suggest that DIC from the active Haima cold seeps probably represent an important source of fossil carbon to the overlying water column. However, more work is required to understand the influence of methane seepage on the chemistry of bottom-water carbonate system, especially the amount of DIC released to the hydrosphere.

Conclusions
The geochemical composition of porewater in shallow sediments and bottom seawater were investigated at four cold seeps in the northwestern SCS. As a result of the study, we conclude the following: 1. Thermogenic or microbial gas transport from deep-seated sediments serves as the main methane source and dissolved at high-flux sites near the seepage centers. Most dissolved methane is consumed by AOM in the cores W08B, W08C, W09, R7-1, R7-3 and HM-1 as indicated by the shallow SMTZ (~7 mbsf to <0.5 mbsf) and our model results. Depth-integrated AOM rates range from 59.0 to 591 mmol m −2 yr −1 . The methane-related carbon turnovers are highly heterogeneous at the studied cold seep systems. We attribute this heterogeneity at cold seeps to the difference in proximity of the sampling sites to the center of fluid conduit.
2. The DIC effluxes range from 51.5 to 427 mmol m −2 yr −1 at the study sites. 13 C-and 14 C-depleted fossil DIC in seep fluids may be released to the bottom water at the four seep areas based on the lower δ 13 C and ∆ 14 C values of bottom-water DIC. Simplified estimations show that cold seeps fluids may contribute 16-26% of the bottom seawater DIC. In addition, a rough average 8.4 × 10 −5 Tmol DIC may be released from shallow sediments to the water column annually at the Haima cold seeps.
Overall, this study shows that the contribution of cold seep fluids is significant both to the pore fluids and the bottom seawater, and may have considerable impacts on the carbon cycle as well as on the seafloor chemosynthetic ecosystems in the cold seep areas.
Since the porosity of our studies cores is not available, we took the mean porosity from an adjacent core [74] and applied it in the model run. As a result, sediment compaction was neglected, and ν s and ν p were equivalent to the sedimentation rate (ω) which is approximated according to [74,75].
Depth-dependent molecular diffusion coefficients of dissolved species were calculated after [73] and [76] and D s were corrected for tortuosity: where D m is the molecular diffusion coefficient in free seawater at the in-situ temperature, salinity and pressure values. The diffusive transport of DIC was simulated using the diffusion coefficient of bicarbonate (HCO 3 − ) since this is the dominant anion.
The major reactions considered in the model are particulate organic matter (POM) degradation via sulfate reduction, methanogenesis, AOM and authigenic carbonate precipitation. The net reaction terms of the one solid (POC) and four dissolved species (SO 4 2− , CH 4 , DIC and Ca 2+ ) are shown in Table A1. All the model parameters are given in Table A2. The length of the simulated cores 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. A zero concentration gradient (Neumann-type boundary) was imposed at the lower boundary for all species. The model was solved using the NDSolve object of MATHEMATICA V. 10.0. All simulations were run for 10 7 years to achieve the steady state with a mass conservation of >99 %. Table A1. Reaction terms of all species used in the model.

Appendix B
Appendix B includes full data of analytical results of the porewater samples of the studied cores (Tables A3 and A4) and the bottom seawater samples (Table A5).