Discharge of deeply rooted fluids from submarine mud volcanism in the Taiwan accretionary prism

Qualitative and quantitative assessments of fluid cycling are essential to address the role and transport of deeply sourced fluids in subduction systems. In this study, sediment cores distributed across a submarine mud volcano (SMV) offshore southwestern Taiwan were investigated to determine the characteristics of fluids generated through the convergence between the Eurasian and Phillippine Sea Plates. The low dissolved chloride concentration combined with the enrichment of 18O, and depletion of 2H of pore fluids suggest the discharge of deep freshwater formed by smectite dehydration at an equilibrium temperature of 100 to 150 °C. The upward fluid velocities, decreasing from 2.0 to 5.0 cm yr−1 at the center to a negligible value at margin sites, varied with the rate and efficiency of anaerobic methanotrophy, demonstrating the impact of fluid migration on biogeochemical processes and carbon cycling. By extrapolating the velocity pattern, the flux of fluids exported from 13 SMVs into seawater amounted up to 1.3–2.5 × 107 kg yr−1, a quantity accounting for 1.1–28.6% of the smectite-bound water originally stored in the incoming sediments. Our results imply that SMVs could act as a conduit to channel the fluids produced from great depth/temperature into seafloor environments in a subduction system of the western Pacific Ocean.

unconsolidated shale/mudstone, fluid overpressurization, and formation of dense fracture array 9,[12][13][14] . Because the fracture network associated with the plate convergence could be extended to as deep as the lower crust, SMVs are able to efficiently export deeply sourced fluids, sediments, and reducing volatiles to the seafloor, generating morphological domes and depressions, mudflows, and even biological hotspots with colonies depending on leaked gaseous hydrocarbons [8][9][10] . While SMVs provide a window to witness the deep subsurface characteristics, upward migrating hydrocarbon gases produced by thermal maturation or methanogenesis also drive the sulfate-dependent anaerobic oxidation of methane (AOM) at shallow depths 9,10 . Such a biological removal mechanism has been considered to be effective in maintaining a low level of methane in seawater 15 . Overall, the exact characteristics and quantities of deep fluids discharged from SMVs, and their impacts on biogeochemical activities remain poorly constrained in the subduction system of the western Pacific Ocean.
In this study, we present a comprehensive geochemical dataset of fluids extracted from sediment cores recovered from sites distributed across an SMV, TY1, in the accretionary prism offshore southwestern Taiwan (Fig. 1). These data were used to assess the fluid source, formation temperature, and fluid-rock ratio for the reaction in the source region. Reactive transport modeling was further applied to quantify the fluxes of fluid and methane exported from various sediment compartments to seawater. Finally, the possible fluid transport pathways in the Taiwan accretionary prism were discussed.

Geological Settings
Offshore southwestern Taiwan is a westward extending accretionary prism related to the subduction of the Eurasian Plate underneath the Philippine Sea Plate since 18 Ma [16][17][18] . The deformation front and its southward extension, the Manila trench, defines the boundary between passive and active margins (Fig. 1a). The active margin is further divided into two structural domains separated by the out-of-sequence thrust: the Lower Slope and the Upper Slope domains 19 (Fig. 1a). The Lower Slope domain is featured with the fold-and-thrust ridges; in contrast, the Upper Slope domain is characterized by mud diapiric structures with mud volcanoes 20 . To date, 13 submarine mud volcanoes have been identified 21 . The largest one, Tsangyao Mud Volcano Group (TYMV; re-named after MV12 22 ), is composed of two mud volcanoes (TY1 and TY2; Supplementary Fig. S1) fed by one mud diapir, as evidenced by seismic reflection profiles 21,23 . The TY1 has a conical structure and a wide flat top with a diameter of ~500 m at a water depth of ~370 m 21 (Fig. 1b,c). Fluid/mud discharges occur at a frequency of every 3 to 10 seconds 24 . Two major gas plumes were detected on the crest by a multibeam echo sounder. These gas plumes could reach to a height of up to 367 m above the seafloor 24 . Other gas venting signals were also observed at the flank of the TY1 ( Supplementary Information, Fig. S2).

Results
Fluid geochemistry. Variation of solute profiles across TY1. Profiles of solute and gas concentrations are shown in Figs. 2 and 3. Most porewater in the upper 120 to 280 cm of sediment had a seawater-like composition. Below this interval at the center sites, concentrations of most ions decreased with depth (potassium, sodium, magnesium, calcium, chloride, and sulfate). In contrast, concentrations of total alkalinity (TA), boron, and lithium increased with depth. Between 120 and 280 cm below seafloor (cmbsf), sulfate was undetectable and accompanied by an increase in methane concentration (1-2 mM), forming a sulfate to methane transition zone (SMTZ). Changes in concentration gradient of other solutes were also observed across the SMTZ. For example, calcium and magnesium concentrations decreased to ca. 1 mM and lower than 5 mM below the SMTZ, respectively. The decrease was coincident with the increase in TA concentration (to around 35 mM).
At the upper flank (site MV12-A) and lower flank (site F6-3) sites, significant decreases in sulfate, calcium, and magnesium concentrations coincided with the increase in methane and TA concentrations (Figs. 2 and 3). Both lithium and boron concentrations remained at constant values throughout the core at the upper flank site. In contrast, boron and lithium concentrations increased with depth at the lower flank site, where SMTZ was observed at the bottom of the core.
At margin sites (sites MV12-C and C-2), concentrations of most ions were close to seawater values. In addition, methane concentrations only increased to at most 50 μM. Both boron and lithium concentrations at site MV12-C decreased with depth (Fig. 3b). , sodium (in green), and potassium (in gray) concentration depth profiles at ten sites investigated in this and previous studies (for sites MV12-A) 22 . The solid lines represent the best fit of the modeled results. The lower boundary conditions at site MV12-1 (in dashed line) were assumed to be the same as those at site A2-2. The modeled upward velocity was also provided. (b) Modeled results for sulfate (blue solid line) and methane (red dashed line) depth profiles and AOM rates (R AOM ; gray area) at eight sites. Hydrocarbon gases and carbon isotopic compositions. At the western and southern center sites, δ 13 C-CH 4 values increased from −40‰ on the seafloor to −20‰ above the SMTZ. The δ 13 C-CH 4 value was the smallest at the SMTZ (−50 to −45‰) and remained at a constant value below the sulfate reduction zone (Fig. 5). At the upper flank (site MV12-A) and lower flank (site F6-3) sites, the δ 13 C-CH 4 value was the smallest at the SMTZ (−70‰ at site MV12-A and −48‰ at site F6-3). The δ 13 C-DIC (dissolved inorganic carbon) patterns for center and flank sites were similar in a way that the values were the smallest at the SMTZ (−30‰ at site 24; −20‰ at sites MD4-P3 and MD4-P1) and increased towards the seafloor (~0‰) and the bottom of the cores.
The C 1 /C 2+ ratios and the δ 13 C-CH 4 values for sites A2-2 and MD4-P3 were between 25 and 90, and between −35.9 and −44.5‰, respectively, falling in the region of thermogenic methane in the plot of C 1 /C 2+ ratio versus δ 13 C-CH 4 value (Fig. 6). For sites 24 and 24-2, the δ 13 C-CH 4 values were between −48 and −37‰, whereas the C 1 /C 2+ ratios increased to more than 120. These data points fell in a region between microbial and thermogenic methane.
Reactive transport modeling and fluid flux from SMVs. The reactive transport modeling yielded that upward fluid velocities at western and southern center sites were in the range of 2 to 3 and 2 to 5 cm yr −1 , respectively (Fig. 2). Upward fluid migration is likely absent at the margin site as the solute concentration do not change significantly with depth. The total fluid discharges from the craters and the rest of the cone structures of 13 SMVs were calculated to be 2.9 to 7.3 × 10 6 kg yr −1 and 1.0 to 1.7 × 10 7 kg yr −1 , respectively (Supplementary Table S6).
The modeled methane fluxes beneath the depth of SMTZ at the western and southern centers were in the range of 1020-1530 and 972-2990 mmol m −2 yr −1 , respectively ( Table 1). The methane fluxes above the SMTZ at these center sites were 162-1380 mmol m −2 yr −1 . The methane flux beneath and above the depth of SMTZ was much lower (347-456 and 9-47 mmol m −2 yr −1 , respectively) at the flank sites (sites MV12-A and F6-3) than at Areas denote with "LRW" and "GW" in (e) marked the ranges of isotopic compositions for local meteoric water and groundwater from adjacent onshore areas, respectively 75 ; gray square area represents the δ 18 O and δ 2 H values of marine clay [48][49][50][51] . Arrows marked with numbers in (d) and (e) represent processes potentially occurring at the source depth, and their directions denote the trends associated with the processes 8 : 1: volcanic ash alteration at temperatures lower than 300 o C; 2: volcanic ash alteration at temperatures higher than 300 o C; 3: gas hydrate dissociation; 4: biogenic opal recrystallization; 5: clay mineral dehydration; 6: meteoric water input. the center sites. While the methane fluxes above and beneath the SMTZ decreased with the distance from the center, the ratios (flux beneath SMTZ/flux above SMTZ) were relatively higher at the flank sites than at the center sites. Such a variation was translated into an increase in AOM efficiency from 54% at the center sites to 98% at the lower flank sites (Table 1).

Discussion
The decreasing chloride concentration with depth ( Fig. 4a) and carbon isotopic composition and abundance patterns of methane (Figs. 5 and 6) suggest an input of freshwater and thermogenic methane from great depth. Crossplots among δ 18 O values, δ 2 H values, and chloride concentrations (Fig. 4e) indicate the deviation of the fluid compositions for the center sites from the GMWL and seawater. The geochemical variation could be best accounted for by the impact of seawater and a deep fluid characterized by low chlorinity, high δ 18 O values, and low δ 2 H values. Such variations in δ 18 O and δ 2 H values cannot be explained by gas hydrate dissociation, ash alteration at high or low temperature, or biogenic opal recrystallization, as these reactions enrich fluids with 2 H or cause no fractionation on hydrogen isotopes 8,10 . Clay dehydration appears to be the most plausible mechanism to account for the observed geochemical signature.
To approximate the temperature of fluid formation, specific cation-ratios were used to derive the fluid temperature [25][26][27] . The validity of geothermometers lies in the requirement of the presence of reactant mineral. Since sediments or mineralogy in the source region are not available, the constituent mineral for the formation currently outcropping in southwestern Taiwan was first adopted to evaluate the applicability of the specific geothermometer. The inference is based on the fact that the subduction of the Eurasian Plate underneath the Philippine Sea   www.nature.com/scientificreports www.nature.com/scientificreports/ Plate generated a series of westward propagating thrusts that sequentially exhumed the sediments once deposited offshore Taiwan 18 . Therefore, the strata on land could mirror the characteristics and constituents of sediments in the potential source region that is accessible only with the deep drilling. The Plio-Pleistocene Gutingkeng Formation formed in shelf break environments was chosen as the analogy for the source material. The formation is primarily composed of mudstone with lens of siltstone and sandstone 28 . Feldspars at abundances of up to 5% in this formation have been observed 29 . In addition to the constraint from the formation in terrestrial environments, feldspars in cores retrieved from the ODP 1144 site in the abyssal region of the South China Sea amounted up to 1.33 wt.% of the total sediments 30 . Therefore, both observations support the possibility that feldspars are available in the source region and serve the source of sodium and potassium in fluids. The geothermometer based on the Na/Li and Na/K ratio were used for the temperature estimation 25 where the concentration is in molality in Eq. (1) and mg/kg in Eq. (2). The Ca-Mg-K and Mg/Li geothermometers were not considered because the porewater chemistry for all depths at the TY1 is oversaturated with respect to a number of carbonate minerals, such as dolomite and Mg-calcite (Supplementary Table S1). Furthermore, the concentration of magnesium is lower than what is expected by the dilution of seawater by freshwater alone, again suggesting the potential removal of magnesium by carbonate precipitation (Supplementary Fig. S3). The SiO 2 geothermometers were not applied due to the potential underestimate caused by rapid cooling through the intrusion of the overlying seawater 26 .
The fluid temperature was calculated to range between 100 and 150 °C. Our estimated temperature is supported by the high concentrations of lithium and boron (Fig. 3b), a pattern consistent with the results from laboratory experiments by which lithium and boron are preferentially released from silicate minerals in marine sediments starting at ca. 50 °C [32][33][34] . The estimated temperature is also consistent with the predominance of thermogenic methane over microbial methane as being inferred on the basis of isotopic compositions and abundance ratios (Fig. 6). Furthermore, the geothermal gradients measured at TY1 vary from 25 °C/km at the margin sites to 390 °C/km at the center sites 23 . Since the crater is directly connected with the conduit channeling hot fluids from great depth, the measured geothermal gradient would be distorted toward a higher value by the advective fluid flow 23,35 . Instead, the averaged value 35 (25 °C/km) measured at the margin sites and in surface sediments in the region was chosen to represent the regional geothermal gradient. Using this value, the depth producing the observed geochemical signatures was calculated to range between 3.6 and 5.7 kmbsf (kilometer below seafloor).
These lines of evidence indicate that a significant fraction of pore fluid is contributed from clay dehydration at high temperatures. Illite and chlorite, which are considerably abundant in sediments offshore southwestern Taiwan 36 , are excluded because illite and chlorite are dehydrated at a much higher temperature range (>237 °C for illite 37 and >650 °C for chlorite 38 ). Kaolinite is also excluded as a result of being less abundant in the northern South China Sea 39 and stable over the temperature range estimated in this study (decomposed at >800 °C and 19 GPa 40 ). Since the transformation of smectite to illite proceeds at a temperature higher than 60 °C 8,34 and produces fluids enriched in 18 O and depleted in 2 H, smectite is considered to be the most plausible mineral candidate to account for the observed geochemical signature. As abundant smectite and mixed layer (1.1-24.8 wt.% in bulk sediments) have been found in the active margin 41 and passive margin off southwestern Taiwan 39,42 , and sediments off the Luzon Arc 43-46 (ODP, sites 1144 and 1146), smectite sourced from Taiwan orogeny and Luzon Arc volcanism and later being subducted/accreted as part of the Taiwan accretionary prism is the most likely source.
To further assess the smectite dehydration, we adopted an isotopic mass balance to calculate the ratio of water to rock (W/R) for reactions considering the exchange of oxygen and hydrogen between minerals and pore fluid under a closed-system condition 47 . The fluid-to-mineral ratio for the reaction was calculated based on the mass balance for hydrogen and oxygen isotope compositions 47 : values were assumed to range from +19 to +28‰ and from −80 to −36‰, respectively. The temperature dependence for the fractionation factors for oxygen and hydrogen isotopes was calculated according to previous studies for a temperature range between 100 and 150 °C 52,53 . The calculations yielded the δ 18 O and δ 2 H values varying as a function of equilibrium temperature and W/R ratio (Fig. 7). Using the measured isotopic compositions of the porewater collected from the center sites (site MD4-P3, Fig. 1b), the W/R ratios were estimated to be <0.5, 0.3-0.9, and 1.1-1.9 if the δ O r i 18 values were set as +19‰, +20 to +22‰, and +23 to +28‰, respectively. The amount of needed smectite could be further deduced from these W/R ratios considering the possible volume of porewater and assuming the complete reactivity of all available porewater and smectite. To constrain the possible porosity at the depth estimated for smectite dehydration, the porosity variation for cores collected from the Nankai Trough 54 (IODP 358, site C0002, 3.2-km long) was used for the extrapolation. The measurement from IODP 358 yielded a porosity of 17% at the depth of 3.2 kmbsf 54 . By applying the exponential reduction function, the porosity at 5.7 kmbsf was calculated to be 5.6% (Supplementary Table S2). Assuming that the porespace is saturated with porewater and that the porosity is 6%, the quantity of porewater would be lower than 3 wt.% of the bulk sediment. To fulfill the estimated W/R ratios (0.3-1.9) for smectite dehydration, 1.6-10 wt.% of smectite is needed. Such a mineral abundance is consistent with the measured amount of smectite in the current active and passive margin 41,43-45 (1.1 to 6.2 wt.% in sediments collected by a sediment trap deployed southwestern Taiwan; 4.5 to 24.8 wt.% in sediments from the ODP 1144 and 1146).
As subducted-smectite could be the source of freshwater, fluid fluxes in the SMV region of the accretionary wedge were further assessed. The total amount of smectite-bound fluid within the incoming plate was estimated considering the weight of the incoming plate and the content of smectite and porosity in the sediments [55][56][57] : where F a is the flux, C a is the weight percentage of mineral-bound water in smectite, H is the pre-subducted sediment thickness, ρ is the density of dry sediments, ϕ a is the porosity of sediments from the incoming plate, L a is the length of trench, and V a is the subduction rate. The details for calculations were stated in Supplementary Information. The total smectite-bound fluid from the incoming plate was calculated to be 0.7-8.8 × 10 8 kg yr −1 .
As sediments are subducted, the buried seawater-like pore fluid and produced smectite-bound fluid would be squeezed and lost gradually with the increasing temperature and pressure. As a result, the actual amount of water released into SMVs region would be substantially lower than that from the incoming plate. Of the overall flux of brackish fluid expelled from SMVs offshore Taiwan (1.3-2.5 × 10 7 kg yr −1 ), approximately 80% of this flux could be attributed to the smectite-bound fluid (1-2 × 10 7 kg yr −1 ) considering the mixing of seawater and freshwater end components (based on the chloride concentration, 110 mM, at the bottom of site MD4-P3). This quantity of smectite-bound fluid released from the SMVs in the region constituted about 1.1-28.6% of the smectite-bound fluid originally subducted into the trench (0.7-8.8 × 10 8 kg yr −1 ).
Previous modeling study using the data integrated from geophysical and mineralogical observations for the Nankai Trough 6 demonstrated that the majority of the smectite-bound fluids (4.17 × 10 8 kg yr −1 ; 88% of total smectite-bound fluids in the incoming plate) would be lost to the frontal prism, leaving a total of 0.56 × 10 8 kg yr −1 smectite-bound fluids further released into the overlying wedge as the descending sediments reach the region beneath the Kumano basin where mud volcanoes are abundant. If such a fluid budget for the Nankai system is directly applied to the Taiwan subduction system, a total of 0.8-10.4 × 10 7 kg yr −1 ( = 0.7-8.8 × 10 8 / (4.17 × 10 8 / 88%/ 0.56 × 10 8 )) smectite-bound fluids would be released from the descending sediments into the wedge in the Upper Slope. This estimated flux is comparable with the flux of freshwater expelled into seawater through SMVs offshore Taiwan obtained in this study (1-2 × 10 7 kg yr −1 ). Whether such an analogy of the fluid budget could be drawn directly remains to be explored even though the plate geometry, heat flux, mineral type and content, convergence rate, and the distance between trench to mud volcanoes are comparable between the Nankai and Taiwan subduction systems (Supplementary Table S3).
The impact of the deep fluid on shallow subsurface microbial processes is reflected in geochemical profiles. The decreasing rate of upward fluid flow from the center towards the margin sites (Fig. 2) varied concomitantly in a similar and an opposite trend with the modeled AOM rates and efficiency, respectively (Table 1 and  Supplementary Table S5). For example, at site MD4-P3 where the upward velocity was the greatest, the AOM rate was the highest among eight sites (1610 mmol m −2 yr −1 ), and the AOM efficiency was the lowest (54%) ( Table 1). In contrast, the upward fluid flow was absent at the lower flank site, whereas the AOM efficiency reached 98% (Table 1). Such an observation suggests that the upward velocity of the fluid is the primary control of AOM rates. As the capacity of AOM is limited, the efficiency of methane removal could not accommodate the extremely high supply of methane at center sites. This pattern is consistent with those obtained in onshore mud volcanoes in Taiwan 22,58 and SMVs elsewhere (e.g., the Hakon Mosby SMV 59 ).

Conclusions
Porewater and gases for cores distributed along a transect of TY1 were systematically analyzed to assess the spatial variation of fluid velocities and biogeochemical activities and further constrain the fluid budget. The co-variation of the isotopic composition and chloride concentration demonstrates that the majority of fluid retrieved from the crater is derived from smectite dehydration at a temperature of 100 to 150 °C or a depth of 3.6 to 5.7 kmbsf. Assuming that the isotopic equilibrium is reached, the water-rock interaction proceeds at a water-rock ratio of more than 0.3. Reactive transport modeling for compounds inert to microbial processes yielded that the velocity of upward migrating fluid varied from 2 to 5 cm yr −1 at the center sites to a negligible level at the margin sites. By extrapolating such a variation pattern for other SMVs offshore Taiwan, the overall flux discharged from SMVs ranged from 1.3 to 2.5 × 10 7 kg yr −1 . This fluid quantity accounts for 1.1-28.6% of the smectite-bound water originally stored in the incoming sediments, implying that SMVs could act as a conduit to channel the fluid produced from great depth/temperature into seafloor environments. The flux variations across the transect of a mud volcano strongly impact the magnitude and efficiency of the biological filtration of methane through methanotrophy. www.nature.com/scientificreports www.nature.com/scientificreports/

Methods and Materials
Sampling sites. Four expeditions were conducted in the offshore southwestern Taiwan by R/V Ocean Researcher 5 (legs 1309-2 and 0039) and R/V Ocean Researcher I (legs 1107 and 1118). Piston and gravity cores with lengths ranging from 102 to 720 cm were recovered along a transect from the center to the margin of the TY1 cone structure (Fig. 1c). Sampling sites are listed in Table 2 and shown in Fig. 1b,c. Sampling and analytical methods. Sediments and core top seawater for gas analyses were collected and transferred into serum bottles. These bottles were later filled with saturated NaCl 60 (in cruises OR5-1309-2 and OR5-0039) or 1 N NaOH 58 (in cruises OR1-1107 and OR1-1118) as preservatives, sealed with butyl rubber stoppers, crimped with aluminum caps, and stored upside down at room temperature or 4 °C. Pore fluid samples for aqueous geochemistry were obtained through centrifugation and subsequently filtered through 0.22-µm nylon membrane syringe filters. The filtrates were further split into five fractions for the analyses of anion, cation, total alkalinity (TA), carbon isotopic compositions of DIC, and oxygen and hydrogen isotopic composition, and stored at 4 °C. For the cation samples, concentrated nitric acid (70%) was added at a volume ratio of 1:45 to preserve the valence state for elements sensitive to the redox change.
The concentrations of headspace hydrocarbon gases were determined by a gas chromatography (GC; SRI 8610 C and Agilent 6890 N) with the precision typically better than 5% 22,58 . Ion concentrations were determined by an ion chromatography (IC, 882 Compact IC Plus) equipped with a Metrosep A Supp 5 column for anions and a Metrosep C 4 column for cations and with the detection limit around 0.1 ppm in weight 61 . Lithium and boron were analyzed by a sector field inductively coupled plasma mass spectrometry (SF-ICP-MS) for 7 Li and 11 B at low resolution (M:ΔM ~300). Six standards were prepared from reference material NASS-5 (all from National Research Council Canada) with the uncertainty better than 3%. TA was determined by titrating the sample with 0.02 N HCl while bubbling with nitrogen 62 . Carbon isotopic compositions of methane and DIC were measured by an isotope ratio mass spectrometry (IRMS) in line with a gas chromatography and a combustion oven. Hydrogen and oxygen isotope compositions of porewater were measured by a LGR Liquid-Water Isotope Analyzer (LWIA; for cruises OR5-1309-2 and OR5-0039) and a Picarro L2140-i Analyzer (for cruises OR1-1107 and OR1-1118). The obtained isotopic compositions were expressed as the δ notation referenced to standards [δ = (R sample /R std − 1) × 1000‰, where R is the ratio of heavy over light isotopes]. Standards are the Vienna Pee Dee Belemnite (VPDB) for carbon isotope, and the Vienna Standard Modern Ocean Water (VSMOW) for hydrogen and oxygen isotopes. The precisions are ±0.3‰ for δ 13 C values by the IRMS measurements, ±0.2‰ for δ 18 O values and ±2.0‰ for δ 2 Η values by the LWIA measurements, and ±0.03‰ for δ 18 O values and ±0.22‰ for δ 2 Η values by the Picarro spectroscopic measurements. The water content of cored sediments was calculated based on the weight loss of the sediments after freeze-drying. The porosity was converted from the water content assuming a dry sediment density of 2.7 g cm −3 63 .

Reactive transport modeling.
To deduce the velocity of upward fluid in such a dynamic environment 64 , a simplified one-dimensional reactive transport modeling was applied for porewater geochemistry profiles with a revised code 65 . The concentration-depth profiles of five dissolved species, including chloride, sodium, potassium, sulfate, and methane, were calculated. This model is formulated by a partial differential equation for solute transport and reactions 65,66 . The modeling was first performed on conservative solutes (chloride and sodium) exempted from microbial processes and most abiotic reactions in shallow sediment environments to determine the range of upward advection rates. The sodium and chloride concentrations in the shallow sediments (the top 120 to 280 cm of sediments; Fig. 2) did not vary much and cannot be attributed to bioirrigation 67 . Therefore, irrigation induced by bubble transport was incorporated to explain such constant concentrations 60,68 . Microbial reactions involving sulfate and methane were considered 12,65 . A detailed description of model-construction, chosen parameters, and boundary condition in the model along with additional sensitivity tests are given in the Supplementary Information. Fluid discharge from SMVs. To estimate the total discharge of brackish water from SMVs, the upward fluid velocities across the TY1 were compiled first with those obtained in the previous study 22 . The upward fluid velocities generally decreased with increasing distance from the edge of crater. Such a trend is also found in other www.nature.com/scientificreports www.nature.com/scientificreports/ SMVs in offshore Norway and the Gulf of Cadiz 12,69 (such as the Håkon Mosby and Carlos Ribeiro SMVs). To generalize the fluid velocities across the entire cone structure, the variations in velocity with distance from the edge of the crater were fitted with the exponential equation ( Supplementary Information, Fig. S9a). The calculation was based on the assumption that the fluid flow was confined through porespace and driven by a single, pressurized fluid source.
The velocity distribution was further integrated with the ring area coverage ( Supplementary Information,  Fig. S9b) to derive the total brackish water discharge using the following equation 12  In this calculation, the crater radius of 250 m, the cone diameter of 3600 m 21 , and the porosity of 36% (an average value of the porosity of TY1; Supplementary Fig. S6) were used. The density of water was assumed to be 1000 kg m −3 . The same approach was applied to other 12 SMVs assuming that the upward velocity at the crater of individual SMVs was the same as that at TY1. Since each SMV has different diameters of its crater and cone structure, the relationship between the velocity and the distance ratio was constructed (normalized to the distance from the edge of crater). Such a relationship was then applied to calculate the discharge fluxes of other individual SMVs ( Supplementary Information, Table S6). Terrestrial mud volcanoes were not considered because their sizes are too small when compared to SMVs 70 . More details are described in Supplementary Information.