MV#1 is in the northern area of the mud volcano field off Tanegashima Island (30° 52.7′N, 131° 46.0′E, 1,400 m water depth; Figures 1A, B); it is ∼2.5 km in diameter and ∼280 m height. During cruise NT15-15 of R/V Natsushima in 2015, we surveyed the summit of MV#1 during two dives (#1866 and #1867) of ROV HYPER-DOLPHIN (Figure 1B). We observed a small colony of Calyptogena (Figures 2A, B). Around the summit, we observed several consolidated sandy mudstone clasts (several to several tens of centimeters in size; Figures 2C, D) that were probably mixed with sediments fractured by diapiric mud intrusion. The presence of these clasts suggests that sediments around the summit comprise erupted mud and not hemipelagic sediments. During cruise KH-15-02 of R/V Hakuho Maru in 2015, two sediment cores were retrieved from MV#1 using a piston coring system: one was used for pore water analysis; the other for sediment analysis.

MV#2 is near MV#1 (30° 55.1′N, 131° 50.6′E, 1,430 m water depth; Figures 1A, B); it is ∼1.2–1.5 km in diameter and ∼120–170 m height. At MV#2, two piston cores and two multiple cores were retrieved for pore water and sediment analyses during cruise KH-19-5 of R/V Hakuho Maru in 2019. Since only conventional piston coring was conducted, a visual survey for MV#2 could not be performed.

MV#3 is in the northern area of the mud volcano field (31° 02.5′N, 131° 41.0′E, 1,200 m water depth; Figures 1A, B); it is ∼1.5 km in diameter and ∼100 m height, making it the smallest of the studied mud volcanoes. During cruise KH-19-5, two piston cores and two multiple cores were retrieved. One of the piston cores was retrieved using the Navigable Sampling System (NSS), which allows observation of the target seafloor area during piston coring via an underwater camera coupled to the piston coring system; the other piston core (for sediment analysis) was retrieved using the normal piston coring system. These camera observations revealed a Calyptogena colony near the summit coring site of MV#3 (Figure 2E), suggesting that this mud volcano was quite active.

MV#14 is situated in the southern area of the mud volcano field (30° 11.5′N, 131° 23.5′E, 1,680 m water depth; Figures 1A, B); it is ∼3 km in diameter and ∼270 m height, making it the largest of the studied mud volcanoes, and it has developed a summit caldera (Figure 1B). During cruise KH-15-2, one NSS piston core was retrieved; in the camera footage, we observed a Calyptogena colony on the eastern edge of the caldera (Figure 2F). Next to the colony, we observed a white, chimney-like object; it may have been a carbonate chimney, but the low resolution of the image precludes a precise interpretation (Figure 2F).

A multiple core was retrieved from a reference site near the mud volcano field during cruise KS-21-18 of R/V Shinsei-Maru in 2021 (31° 29′N; 132° 02′E, 1,220 m water depth; Figure 1A).

Piston core samples retrieved for sediment analyses [Length: 3.11–5.21 m, inner diameter: 75 cm (Supplementary Table S1)] were cut into 1-m-long sections and stored onboard at 4°C. Onshore laboratory, the sections were cut in half lengthwise and the sediments were subsampled from the cut surface after the scan for X-ray computed tomography. For the MV#14 core, sediments for pore water and microbial analyses were collected from the piston core onboard, and sediment analysis samples were later subsampled from the microbial analysis samples in onshore laboratory. For nannofossil observations and beryllium isotopic (¹⁰Be) and mineral compositional analyses, sediments were subsampled from the muddy matrix representing the mud volcano sediments; we avoided mud breccias because they were probably incorporated during mud diapir intrusion (Ujiié, 2000).

To observe the internal structure of the core sample, we obtained X-ray computed tomography (X-ray CT) image using Aquilion PRIME/Focus (Canon Medical Systems Corporation, Japan). For the MV#3 core, the most sediments were collected for geochemical and microbial analyses onboard. Therefore, the residual sediments on the core-liner were scanned. For the MV#14 core, all sediments were sub-sampled on board, so we were unable to scan the core.

To examine the origin of the source sequence of the mud volcano sediments, the age of the mud volcano sediments was estimated based on the calcareous nannofossils in the sediments. The samples were prepared as smear slides using standard techniques (Bown and Young, 1998), and then examined with a Nikon E600 polarizing light microscope at 1,500× magnification. Species identification and their occurrence age range of each species followed Backman et al. (2012) and Nannotax3 (https://www.mikrotax.org/Nannotax3/index.html). Semi-quantitative data were collected from the samples, which contained sufficient number of calcareous nannofossils, by identifying at least 200 upper photic zone specimens in a varying number of fields of view per sample. Lower photic zone species Florisphaera profunda was counted separately in the same fields of view when encountered. Upper photic-zone species were classified into six groups (Emiliania huxleyi, Pliocene–Pleistocene taxa except E. huxleyi, Miocene–Pleistocene taxa, Paleogene–Neogene taxa, Miocene taxa, and miscellaneous) based on their occurrence age range (Supplementary Table S2), and the relative abundance of each group was calculated (Supplementary Table S3).

The age of the source sequence of the mud volcano sediments was also estimated using the beryllium isotope analysis. Meteoric ¹⁰Be is formed in the atmosphere through cosmic ray interactions with oxygen and nitrogen nuclei, eventually reaching the ocean via precipitation and aerosols (e.g., Brown, 1987). In the ocean, it is incorporated into seafloor sediments through scavenging by biogenic opal, clay, carbonate, and organic materials (e.g., McHargue and Damon, 1991). It undergoes radioactive decay with a half-life of 1.36 million years, transforming into stable ¹⁰B [Nishiizumi et al. (2007)]. Therefore, unless sediments are exposed at the surface after the deposition, their age can be can be estimated by the concentration of ¹⁰Be in sediments (Amin et al., 1975).

The 0.3 g of sediment samples were decomposed at 110°C for 12 h in a tightly sealed vessel with 1 mL of HF, 6 mL of HNO₃, and 1 mL of HClO₄. After decomposition, 0.3 mg of Be carrier was added to the sample, which was then dried on a hotplate and converted to Cl manner by repeated drying with 2 mL of 35% HCl. To mask other interfering elements such as aluminum for the accelerator mass spectrometer (AMS) measurement, 25 mL of 21% EDTA solution with pH 6.0 was added to the sample, followed by the addition of 2 mL of acetylacetone to form the beryllium acetylacetonate complex. The Be complex was then extracted into 5 mL of CCl₄ using solvent extraction, and back-extracted into 10 mL of 1 M HCl. The extracted Be solution was purified using cation exchange with a 5 mL Dowex 50 W-X8 column (inner diameter 9 mm, length 10 mm, 1 column volume = 5 mL) with 1.0 M HCl. The Be fraction, eluted between 5 and 9 column volumes, was concentrated to 1 mL and transferred to a 15 mL polycarbonate spitch centrifuge tube. The purified Be solution was then precipitated as Be (OH)₂ by adding an ammonium solution. To remove ¹⁰B, which causes isobaric interference in AMS measurement, the precipitation was washed three times with Milli-Q water (>18 MΩ), and re-precipitated using ultra-pure HCl and ammonium solution. After washing, the solution was dissolved in a few drops of ultra-pure HCl and transferred to a quartz tube for drying. The purified Be (OH)₂ was converted to BeO by igniting at temperatures above 800°C. The resulting BeO was mixed with Nb powder and pressed into an OFCu cathode with an inner diameter of 1 mm. The reactive ¹⁰Be abundances were measured using a National Electrostatic Corporation (NEC) accelerator mass spectrometer (AMS) with a 4.8 MV terminal voltage at the Micro Analysis Laboratory, Tandem accelerator (MALT) of the University of Tokyo (U Tokyo), following the method described by Matsuzaki et al. (2020). Absolute values were obtained using the KNB5-2 standard (¹⁰Be/⁹Be = 8.56 × 10⁻¹²; Nishiizumi et al., 2007) with a typical beam current of 4 micro A (⁹Be¹⁶O; Yokoyama et al., 2019). The ratios of ¹⁰Be/⁹Be in the samples ranged from (7–20) × 10⁻¹⁴, while the carrier blank had a ratio of 1.4 × 10⁻¹⁴.

To reveal the mineralogical characteristics of the mud volcano sediments, whole-rock and clay-fraction mineral compositions were determined by XRD using an X’Pert PRO MPD X-ray diffractometer system (Malvern Panalytical, Ltd., Malvern, United Kingdom) with CuKα radiation at 45 kV and 40 mA, 0.0083° divergence, anti-scattering slits, and a 0.15 mm receiving slit. The diffracted X-rays were counted for 10.2 s at steps of 0.0083° 2θ. For whole-rock and clay-fraction mineral compositions, we scanned 3°–75° 2θ and 3°–35° 2θ, respectively.

The whole rock mineral compositions were quantified based on relative peak intensities of XRD pattern and the Rock jock program (Eberl., 2003). The analysis procedure followed Tanikawa and Tadai (2018). Sediment samples were dried at 50°C, and then powdered using an agate mortar. Then, 1 g of powdered sample was mixed with 0.25 g of high-purity Al₂O₃ (AX-5HM, Hinomoto Kenmazai CO., Ltd., Tokyo, Japan), which is the internal standard for Rock jock analysis, using an agate mortar. The 0.1 g of the mixed sample was placed on a slide glass and subjected to XRD analysis.

For clay-fraction analyses, the clay size fraction (<2 μm) was separated from the bulk sediment samples by centrifugation. Clay suspensions were dropped onto glass slides and dried in an oven at 50°C to prepare oriented samples for XRD analysis. After a first round of XRD analysis, the samples were saturated with ethylene glycol vapor at 50°C for more than 5 h and reanalyzed. We also semiquantitatively estimated the proportion of major clay minerals in the clay size fraction using the method of Biscaye (1965). The proportions of illite and smectite crystallites in illite/smectite (I/S) mixed-layer clays were calculated using the angular separation (Δ° 2θ) between the composite illite-(001)/smectite-(002) reflection at 2θ ≈ 10.5° and the composite illite-(002)/smectite-(003) reflection at 2θ ≈ 16° (Moore and Reynolds, 1997).

The percentage of the clay size fraction in the bulk rock was determined by particle size analysis using a Mastersizer et al., 2020 (Malvern Panalytical, Ltd., Malvern, United Kingdom). The sediment samples were first washed through a 2-mm sieve with Milli-Q water, stirred in a dispersant (0.05% sodium hexametaphosphate solution) for 24 h, and then analyzed.

The vitrinite reflection was analysed to examine the thermal maturity of organic matter. The vitrinite reflectance value (R_O) increases exponentially with temperature and records the highest temperature experienced by organic matter (Burnham and Sweeney, 1989). We collected samples of organic matter from MV#1, #2, #3, and #14. To collect enough samples for the vitrinite analysis, sediment samples were collected from three depths in each core, and the data was combined for each mud volcano. For each sample, coal fragments were separated and concentrated using heavy liquids (a sodium polytungstate solution with a specific gravity of 1.8). The separated grains were then mounted in resin and polished using 0.06-μm alumina powder. The measurements were conducted with a silicone-diode microphotometer using 546 nm non-polarized light, in line with the regulations of the International Organization for Standardization (ISO 7404-5), the American Society for Testing and Materials (ASTMD 2798-99), and the Japan Industrial Standard (JIS M 8816). A concentrated spot beam of 1.6-μm diameter was used to measure the reflectance. We measured mean random reflectance (R_m) following ASTMD 2798-05. A linear function between mean random reflectance and mean maximum reflectance has been reported in case of lower vitrinite reflectance (<3.0%: Chijiwa, 1990; <1.7%: ASTM Standard D2798-05). Therefore, the effect of anisotropy of vitrinite reflectance should be small when it was low. A total of >41 grains per each sample were measured in the immersion oil which has a refractive index of 1.516, where a sufficient number of polished grains were present on the sample surface. The equations for calculating the maximum temperature experienced by the sample using measured R_m values were defined by Sweeney and Burnham (1990).

X-ray computed tomography (CT) images of the cores from MV#1–3 showed that mud-supported breccias are included in the muddy sediments (X-ray intensity being proportional to density; Figure 3). The observed lack of sedimentary structure suggests that clasts of surrounding rocks were ripped away from their source strata and incorporated into the intruding mud diapir via hydraulic fracturing (Ujiié, 2000). Sediments from MV#2 in particular contain relatively hard, centimeter-scale clasts, especially at greater depths. However, these clasts in Section 5, Section 6 are loosely consolidated, so that they crumbled easily in hand. Note that MV#3 core was only imaged for less than half of its width due to scanning residual sediments after on-board sub-sampling. In contrast to the mud volcano sediments, the reference core did not contain any evidence of a mud breccia.

All studied samples from the cores MV#1, MV#2, and MV#14 yielded sufficient number of calcareous nannofossils for sub-quantitative study. Samples from the core MV#3, however, varied in abundance of calcareous nannofossils. The samples from 47.5, 153.5, and 314.5 cm below seafloor (cmbsf) contained calcareous nannofossils abundantly, while the samples from 126.5 and from 263 cmbsf rarely yielded calcareous nannofossils (Supplementary Table S2).

The calcareous nannofossil assemblages of all studied samples contained long-ranged species such as the Eocene–Pleistocene species Reticulofenestra minuta and/or the Miocene–Pleistocene medium-sized Reticulofenestra species (Reticulofenestra minuta, Reticulofenestra minuta/haquii, Reticulofenestra pseudoumbilicus (5–7 µm), or Reticulofenestra producta) common-abundantly, as well as various Paleocene–modern species such as Coccolithus pelagicus in relatively low abundance. Twelve Miocene–specific species were recognized in low abundance in many samples, but the Paleocene–specific species were not recognized (Figure 4; Supplementary Table S2). These results suggest that the age of sediments of source sequence of mud volcanos is in the Miocene. Calcareous nannofossil assemblages from the uppermost sample (22 cmbsf) at MV#14 was dominated by Emiliania huxleyi, and lacked Miocene-specific species, suggesting that MV#14 is currently inactive and the uppermost sample (22 cmbsf) of the site was deposited after the end of mud–eruption activity in/before the E. huxleyi acme period (0.082 Ma–present, Raffi et al., 2006). Coexistence of E. huxleyi and Miocene–Pleistocene species in upper 40 cmbsf at MV#1 and around 120 cmbsf at MV#14 implies the incorporation of the Quaternary sediments by the movement of mud associated with the occasional mud eruption (Fesker et al., 2014). Ujiié (2000) reported that the Quaternary hemipelagic sediments, approximately 3.5 m thick, overlay the mud volcano sediments in MV#2, indicating that mud volcanism might have been inactive at the time of their 1994 survey. However, the absence of E. huxleyi in the surface sediments from MV#2 collected during 2019 survey suggests that mud eruptions may have occurred after 1994 or that the 1994 coring missed the location where mud was erupting.

The contents of ¹⁰Be varied over a wide range across the four studied mud volcanoes (6–318 × 10⁵ atoms g⁻¹; Figure 5A; Supplementary Table S4). At MV#1, ¹⁰Be contents were relatively constant at ∼115 × 10⁵ atoms g⁻¹ regardless of depth. At MV#3, ¹⁰Be contents were ∼25–50 × 10⁵ atoms g⁻¹. At MV#14, ¹⁰Be contents were 144–318 × 10⁵ atoms g⁻¹, higher than those at the other mud volcanoes, and were greatest in the shallowest sample (22 cmbsf). At MV#2, we analyzed ¹⁰Be contents in the mud breccia (>2 mm), mud matrix (2 μm–2 mm), and clay (<2 µm) fractions separately to test for any size-dependence of ¹⁰Be contents (ages). Ujiié (2000) reported that the exotic fossils and clast shapes in the mud volcano sediments were torn off and taken into the mud diapir via hydraulic fracturing through an intrusive process. Therefore, we expected the ¹⁰Be contents are different among the mud breccia, mud matrix, and clay. However, the ¹⁰Be contents were not consistent across all particle sizes at each depth. Breccia, matrix, and clay ¹⁰Be contents at 43 cmbsf were 22 × 10⁵, 6 × 10⁵, and 19 × 10⁵ atoms g⁻¹, respectively (lowest content in the matrix, highest in the breccia); those at 216 cmbsf were 16 × 10⁵, 34 × 10⁵, and 30 × 10⁵ atoms g⁻¹, respectively (lowest in the breccia, highest in the clay); and those at 433 cmbsf were 147 × 10⁵, 272 × 10⁵, and 261 × 10⁵ atoms g⁻¹, respectively (lowest in the breccia and highest in the clay). More data would be required to investigate the ¹⁰Be age differences depend on sediment type.

To calculate ¹⁰Be ages, we used ¹⁰Be data for surface sediments from two sites as initial concentrations: our reference site, off Tanegashima Island (4.4 × 10⁸ atoms g⁻¹), and one in Sagami Bay, south of Honshu, central Japan (14.0 × 10⁸ atoms g⁻¹), because of the wide variation in ¹⁰Be contents in surface sediments ranging from 1.2 × 10⁸ and 16.0 × 10⁸ atoms g⁻¹ (Bourles et al., 1989). The calculated ¹⁰Be ages using the ¹⁰Be content at the reference site ranged from 7 to 15 Ma, and those using the Sagami Bay data ranged from 5 to 13 Ma (Figure 5B). Although the estimated ages deviated by 2 Myr because of the difference in the assumed initial concentration, they are all middle to late Miocene.

This ¹⁰Be age is consistent with the nannofossil assemblages which suggest that the age of the source sequence is the Miocene. Because the nannofossil observations indicate that the sediments comprise a mixture of sediments deposited at different ages, the ¹⁰Be ages should be mixed values and the observed age variations are likely due to differences in the mixing ratio. Indeed, we obtained the youngest ¹⁰Be age (7–5 Ma) for the highest ¹⁰Be content observed at 22 cmbsf in MV#14, consistent with the dominance of E. huxleyi in the nannofossil assemblage at that depth. Furthermore, the estimated ¹⁰Be ages for MV#1 and MV#14 (<10 Ma) were younger than those for MV#2 and MV#3 (>10 Ma), consistent with the higher abundances of Pliocene–Pleistocene taxa except E. huxleyi at MV#1 and MV#14 than at MV#2 and MV#3 (Figure 4). Thus, we attribute the young ¹⁰Be ages at MV#1 and MV#14 to increased mixing with younger sediments. We note that, at 432.5 cmbsf in MV#2, the ¹⁰Be age was relatively young (8–6 Ma). However, Minylitha convallis (Zones NN9–NN11A; ∼9.8–7.7 Ma) was abundant in the same sample (Supplementary Table S2), indicating good agreement between the ¹⁰Be and nannofossil ages. Based on the nannofossil and ¹⁰Be data, we conclude that each studied mud volcano has a middle to late Miocene source sequence. Note that depending on the proportion of younger sediments, there is a possibility that the sediments could attributed to the Early Miocene.

The mud volcano sediments comprised 17–33 wt% quartz (avg. 26 wt%), 14–33 wt% illite (avg., 23 wt%), 12–19 wt% plagioclase (avg., 16 wt%), and other minerals (Figure 6A; Supplementary Table S5). Overall clay mineral contents were 32–52 wt% (illite, smectite, chlorite, and kaolinite; avg., 43 wt%), MV#1 sediments contained less pyroxene (<2 wt%) than other mud volcano sediments (3–7 wt%). MV#1 and MV#14 contained slightly more quartz (25–33 wt%) and slightly less clay minerals (32–47 wt%) than MV#2 and MV#3 (17–24 wt% quartz, 41–52 wt% clay minerals). The generally similar mineral contents of the sediments at each mud volcano (i.e., dominated by clay minerals, especially illite) implies that the mud volcanoes off Tanegashima Island are rooted in the same source sequence.

Figure 7 shows X-ray diffractograms for ethylene-glycolated samples from MV#1–3, MV#14, and the reference site; those of the clay size fraction showed similar patterns. The X-ray diffractograms of the mud volcano sediments showed high peak intensities at 2θ ≈ 9°, which represents illite (001). A small peak at 2θ ≈ 5.5°, attributed to a mixture of discrete smectite and I/S mixed-layer clay, was hidden behind the chlorite (001) peak at 2θ ≈ 6°. The overlapping peaks of chlorite (002) and kaolinite (001) at 2θ ≈ 13.5° were high. In the samples from deeper than 337 cmbsf in MV#2, the peak of discrete smectite and I/S mixed-layer clay was high. In that same sample, we could also identify the composite reflections of illite (001)/smectite (002) at 2θ ≈ 10° and of illite (002)/smectite (003) at 2θ ≈ 16°. These observations indicate that illite dominates the mineralogy of the clay fraction, whereas smectite represents only a small proportion. Peak intensities in the X-ray diffractogram of the reference site sample were much lower than those of the mud volcano sediments (Figure 7); in particular, the smectite peak at 2θ ≈ 5.5° was subtle.

Particle size analysis showed that the clay size fraction (<2 µm) comprises 19%–28% of the mud volcano sediments and 4%–11% of the reference sediment (Supplementary Table S6).

The relative proportions of clay minerals in the clay size fraction, estimated following Biscaye (1965), are shown in Figure 6B and listed in Supplementary Table S6. Below 58 cmbsf at MV#3, the smectite peak could not be detected; therefore, we did not estimate relative clay mineral abundances at greater depths there. In most samples from MV#2 (<337 cmbsf), MV#3, and MV#14, the clay mineralogy of the clay fraction comprised <2% smectite, 65%–79% illite (avg., 74%), and 20%–34% chlorite + kaolinite (avg., 26%). At >337 cmbsf in MV#2, the smectite content increased to as high as 43%. Sediment at the reference site comprised ≲1% smectite, ∼76% illite, and ∼23% chlorite + kaolinite. In MV#3, the integrated peak areas for illite and chlorite + kaolinite peaks below 58 cmbsf are similar to those above 58 cmbsf (Figure 7; Supplementary Table S6). Therefore, it can be inferred that the abundance ratios of illite to chlorite + kaolinite below 58 cmbsf are similar to those above 58 cmbsf.

Overall, the clay fraction at each mud volcano was characterized by similar illite-dominated clay mineralogies. Alongside the similar whole-rock mineral compositions of sediments at each mud volcano, this finding further suggests that the mud volcanoes off Tanegashima Island are rooted in the same source sequence.

Histograms of vitrinite reflectance at each mud volcano peaked at 0%–0.5%, with consistent median values at MV#1, MV#2, MV#3, and MV#14 of 0.42%, 0.43%, 0.44%, and 0.42%, respectively (Figure 8). Considering the thresholds for oil, wet gas, and dry gas production (0.5%–0.6%, >1.3%, and >1.7%, respectively; e.g., Tissot et al., 1987), these low reflectance values indicate that the organic matter transported to the surface at these mud volcanoes is too immature to produce oil or gas. The similar reflectance values and distributions at each mud volcano also indicate that the organic matter is derived from the same source sequence. Assuming the age of the sediment to be 15 Ma (i.e., our oldest ¹⁰Be age), the experienced temperatures of organic matters for MV#1, MV#2, MV#3, and MV#14 were calculated to 81, 83, 85, and 81°C, respectively (Supplementary Table S7).

The calcareous nannofossil assemblages, ¹⁰Be ages, similar illite-dominated mineralogical compositions, and vitrinite reflectance values (∼0.4%) of the sediments at the studied mud volcanoes strongly suggest that the mud volcano field off Tanegashima Island is rooted in the same Miocene source sequence containing immature organic matter. Although sediments at MV#1 and MV#14 are characterized by lower abundances of Miocene nannofossils, younger ¹⁰Be ages, and lower clay mineral contents than MV#2 and MV#3, we attribute these differences to mixing with younger sediments.

Based on the analyzed sediment data, we can constrain the origin of the mud volcano sediments. The Upper Shimanto Group, characterized by Paleogene and lower Miocene rocks, is distributed along Tanegashima Island and the southern parts of Kyushu and Shikoku Islands (Taira et al., 1982). Thus, the group should also occur in the subseafloor in the mud volcano field off Tanegashima Island. If the age of the source sequence is early Miocene, i.e., older than our estimation, the source sediments might be accreted sediments of the Upper Shimanto Group. However, if the source sequence is middle to late Miocene as estimated here, the source sequence cannot be derived from the Shimanto Belt. Based on the experienced temperatures determined from the vitrinite reflectance values (81°C–85°C) and the local geothermal gradient of 25°C km⁻¹ (Ohde et al., 2018), we estimate the depth of the source sequence to be 3.2–3.4 kmbsf.

The illite-rich sediments in the studied mud volcanoes off Tanegashima Island contrast with the smectite-rich sediments in mud volcanoes of the Kumano forearc basin (Ijiri et al., 2018b). The difference could be explained by the distribution of illite in the present surface sediments. Southeast of Kyushu Island, including our study area, illite contents are extremely high in the clay fraction of surface sediments (Aoki and Oinuma, 1988). The illite is detrital in origin, supplied to the East China Sea from mainland China by the Yangtze and Yellow Rivers (Nittrouer et al., 1984; Aoki and Oinuma, 1988) and transported from the East China Sea by the Kuroshio Current (Yin et al., 1987; Aoki and Oinuma, 1988). If the source sequence was deposited during the Miocene, the East China Sea was not yet fully open during their deposition (e.g., Fu et al., 2022), and the illite in those sediments may have been transported directly from mainland China. Nonetheless, this variation in supply paths cannot explain the high clay mineral contents of the studied mud volcano sediments.

Cl⁻ concentrations in pore waters at all studied mud volcanoes decreased from the seawater concentration (∼560 mM) with increasing depth (Figure 9A; Supplementary Table S8). The Cl⁻ concentrations decreased asymptotically to ∼230 mM by 250 cmbsf at MV#1, ∼250 mM by 260 cmbsf at MV#2, and ∼220 mM by 80 cmbsf at MV#3. In contrast, the Cl⁻ concentration only decreased to ∼480 mM at the core bottom (294 cmbsf) at MV#14. The general asymptotic decrease of Cl⁻ concentrations with increasing depth at all studied mud volcanoes suggests the upward advection of low-Cl⁻ fluids.

Porewater δD and δ¹⁸O profiles were consistent with the Cl⁻ profiles (Figures 9A–C; Supplementary Table S8). The δD values decreased asymptotically from ca. −2‰ near the seafloor to ca. −9‰ by 250 cmbsf at MV#1, ca. −15‰ by 260 cmbsf at MV#2, ca. −15‰ by 80 cmbsf at MV#3, and ca. −4‰ to −3‰ by 70 cmbsf at MV#14. The δ¹⁸O values increased asymptotically from ∼0‰ near the seafloor to ca. +5‰ by 250 cmbsf at MV#1, ca. +3.5‰ by 260 cmbsf at MV#2, ca. +3‰ by 80 cmbsf at MV#3, and +1.1‰ to +1.6‰ by 90 cmbsf at MV#14.

Plots of δD and δ¹⁸O values against Cl⁻ concentrations reveal linear relationships at each mud volcano (Figures 9D, E), with the exception of δD vs. Cl⁻ at MV#14 (R ² = 0.01). Enrichment in ¹⁸O and depletion in D alongside depletion in Cl⁻ is attributed to the incorporation of freshwater due to clay-mineral dehydration, such as during the smectite-to-illite transformation that typically occurs at 60°C–160°C (Yeh, 1980; Kastner et al., 1991; You and Gieskes, 2001). Indeed, this process was previously reported at MV#1 by Nakayama et al. (2010).

Assuming that the porewater isotopic compositions result from simple mixing with freshwater [(Cl⁻) = 0 mM], the δ¹⁸O and δD values of the freshwater component can be estimated by extrapolating from [Cl⁻] = 0 mM to the observed porewater data. We obtained extrapolated freshwater δD values at MV#1, #2, and #3 of −14‰ ± 1‰, −24‰ ± 2‰, and −25‰ ± 1‰, respectively (95% confidence interval); we did not extrapolate the δD value for MV#14 because δD was not correlated with Cl⁻. We obtained extrapolated freshwater δ¹⁸O values at MV#1, #2, #3, and #14 of +8.5‰ ± 0.3‰, +6.3‰ ± 0.4‰, +5.3‰ ± 0.4‰, and +8.0‰ ±3.0‰, respectively (95% confidence interval). The δ¹⁸O values of the freshwater components are similar at MV#1 and MV#14 (+8.5‰ and +8.0‰, respectively) and at MV#2 and MV#3 (+6.3‰ and +5.3‰, respectively). These differences in the end member freshwater values may be due to differences in the δ¹⁸O values of the source clays and/or variations in the temperature at which they dehydrated.

To constrain the reason for the different end member freshwater δ¹⁸O values at each mud volcano, we calculated the formation temperature of the fluid following the method of Hensen et al. (2004). Using the estimated δD and δ¹⁸O values of the freshwater component and previously reported clay–water equilibrium fractionation factors for δD (Yeh, 1980; Capuano, 1992) and δ¹⁸O (Sheppard and Gilg, 1996), we estimated theoretical clay mineral δD and δ¹⁸O values capable of producing the estimated δD and δ¹⁸O values of the freshwater component at different temperatures (Figures 10A, B). Then, we compared previously reported H and O isotopic variations in marine clays (shaded ranges in Figure 10; Savin and Epstein, 1970; Suchecki and Land, 1983; Yeh, 1980) to the theoretical δD and δ¹⁸O values calculated for clay minerals at different temperatures to constrain the fluid formation temperature (Figures 10A, B). The equilibrium fractionation curves for MV#1–3 and MV#14 intersect the δ¹⁸O range of marine clays at 90, 70, 60, and 85°C, respectively, representing minimum values (Figure 10A). We could not estimate maximum values based on the δ¹⁸O range because the high-temperature intersections were around 400°C, at which temperature clay-mineral dehydration should be complete. Nonetheless, we constrained the maximum formation temperature at MV#1, #2, and #3 to 145, 185, and 190°C, respectively, using the δD fractionation factor of Capuano (1992) (Figure 10B). However, that fractionation factor gives low-temperature intersections below 60°C, at which temperature clay-mineral dehydration does not occur. Therefore, we estimate that the clay minerals dehydrated at 60°C–190°C, which is generally consistent with the smectite–illite transformation temperature (60°C–160°C; Yeh, 1980; Kastner et al., 1991) and estimated equilibrium temperatures of heavy noble gases at MV#1 and MV#14 (83°C–230°C; Mitsutome et al., 2023). Given the local geothermal gradient of 25°C km⁻¹ (Ohde et al., 2018), we estimate that clay-mineral dehydration occurred at 2.4–7.6 kmbsf.

To quantitatively evaluate the observed differences in the vertical Cl⁻ profiles, we modeled them down to 10,000 cmbsf in 10-cm increments by using the one-dimensional unsteady diffusion/advection tracer equation: ∂ C ∂ t = ∂ ϕ D e f f ∂ C ∂ z ϕ ∂ z   – ∂ ϕ U C ϕ ∂ z (1) where C is the porewater Cl⁻ concentration, ϕ is porosity, D _eff is the effective diffusion coefficient of Cl⁻ in the sediment with respect to tortuosity, and U is the pore fluid velocity. We assumed ϕ = 0.5, the average porosity of MV#1 sediments. We used the relationship (Van Loon and Mibus, 2015) D e f f = D 0 ∙ ϕ m 1 + B ∙ ϕ m 2 (2) where D ₀ = 2 × 10⁻⁹ m² s⁻¹ is the diffusion coefficient of Cl⁻ in water and m ₁ = 2.4, m ₂ = 1, and B = 10⁻¹¹ are constants, to obtain D _eff = 4.3 × 10⁻¹⁰ m² s⁻¹. As the initial state of the Cl⁻ profile just after a mud eruption, we assumed the Cl⁻ concentration at the bottom of the MV#1 and MV#3 cores (∼230 mM) to be the original value of the deep fluid, and that the Cl⁻ concentrations were constant from depth to the surface immediately after the eruption (dashed lines in Figure 11). Then, we estimated fluid advection rates and the time since the most recent mud eruption by comparing the numerically simulated and observed Cl⁻ profiles. We obtained respective advection rates and times since the most recent eruption of 10 mm y⁻¹ and 100–500 years at MV#1, 1 mm y⁻¹ and 100–200 years at MV#2, 100 mm y⁻¹ and 4–20 years at MV#3, and 0.1 mm y⁻¹ and 2,000–10,000 years at MV#4. At MV#3, the simulated Cl⁻ profiles change little from 4 to 20 years, indicating that the fast advection rate there mainly determines the vertical profile at steady state. In contrast, the slow advection rate and downward diffusion of seawater from the seafloor mainly determine the profile at MV#14. These results indicate that the relative levels of activity at each mud volcano are, in order of decreasing activity, MV#3, MV#2, MV#1, and MV#14, which is dormant. This result is consistent with our nannofossil observations and the conclusions of Mitsutome et al. (2023). The Cl⁻profile in MV#1 sediment obtained during the 2002 survey (Nakayama et al., 2010) is similar to that obtained in this study during the 2015 survey, suggesting the activity of MV#1 has remained unchanged from 2002 to 2015.

Porewater SO₄ ²⁻, CH₄, and DIC concentrations and δ¹³C-CH₄ and δ¹³C-DIC values are shown in Figures 12A–D and reported in Supplementary Tables S8, S9. At MV#1, pore water contained <1 mM CH₄ at 40 cmbsf and >5 mM below 90 cmbsf. The δ¹³C-CH₄ values increased from −58‰ at 40 cmbsf to −42‰ at 290 cmbsf. Porewater SO₄ ²⁻ concentrations decreased from 25 mM at 25 cmbsf to 3 mM at 75 cmbsf and <2 mM below 125 cmbsf. Based on the CH₄ and SO₄ ²⁻ profiles, we interpret that the sulfate–methane transition zone (SMTZ), where most methane is consumed by microbial anoxic methane oxidation and microbial sulfate reduction) at MV#1 is at 70–90 cmbsf. DIC concentrations generally increased from 3 mM at 25 cmbsf to ∼45 mM below 270 cmbsf, and δ¹³C-DIC values increased asymptotically from ca. −20‰ near the surface to +1‰–2‰ below 220 cmbsf.

At MV#2, CH₄ concentrations were 5–10 mM and δ¹³C-CH₄ values varied from −56‰ to −52‰, respectively, largely regardless of depth. SO₄ ²⁻ concentrations were ∼28 mM near the seafloor and decreased to <1.5 mM below 60 cmbsf. The CH₄ and SO₄ ²⁻ profiles indicate that the SMTZ is within the upper 60 cmbsf. DIC concentrations varied within 10–18 mM and δ¹³C-DIC values increased asymptotically from 0‰ at 60 cmbsf to +7‰ below 210 cmbsf.

At MV#3, CH₄ concentrations were 3–7 mM and δ¹³C-CH₄ was stable at ca. −52‰ regardless of depth. SO₄ ²⁻ concentrations were 28 mM near the seafloor, 1 mM at 13 cmbsf in the piston core, and 0–2 mM at greater depths, placing the SMTZ near the seafloor. DIC concentrations generally increased from ∼10 mM near the seafloor to 17 mM at the bottom of the core, except for an exceptionally high value of 25 mM at 130 cmbsf, and δ¹³C-DIC values increased asymptotically from +3‰ near the seafloor to ca. +20‰ below 130 cmbsf.

At MV#14, CH₄ concentrations were <1 mM in the upper 70 cmbsf and increased to 2–5 mM below 90 cmbsf (except for a deviation to 0.7 mM at 219 cmbsf). δ¹³C-CH₄ values decreased from −63‰ at 22 cmbsf to ca −76‰ at greater depths (again, except for a deviation to −57‰ at 219 cmbsf). SO₄ ²⁻ concentrations decreased from 28 mM at 20 cmbsf to <3 mM below 95 cmbsf. The CH₄ and SO₄ ²⁻ profiles place the SMTZ at 70–90 cmbsf. DIC concentrations generally increased from 3 mM near the seafloor to ∼18 mM at the bottom of the core. The δ¹³C-DIC values decreased from −5‰ near the seafloor to −34‰ at 94 cmbsf, then increased to −17‰ at the bottom of the core.

Based on the methane/ethane concentration ratios (C₁/C₂) and δ¹³C-CH₄ and δD-CH₄ values below the SMTZ, hydrocarbon gases at MV#1 are mostly derived from the thermal decomposition of organic matter in deep sediments (Figures 13A, B). At MV#14, these data indicate that most methane is of biogenic origin, i.e., produced by hydrogenotrophic methanogenesis. Based on the C₁/C₂ ratios and δ¹³C-CH₄ values at MV#1 and MV#14, we developed a two-end member mixing model between biogenic and thermogenic methane (Figure 13A). As the biogenic end member, we used the sample with the highest C₁/C₂ value (4,414) and the lowest δ¹³C-CH₄ value (−76‰) from MV#14. We used two different thermogenic end members: the lowest C₁/C₂ ratio (8; associated δ¹³C-CH₄ = −47‰) and the highest δ¹³C-CH₄ value (−42‰; associated C₁/C₂ = 20) from MV#1. This mixing scenario constrains the thermogenic contribution to 60%–80% of the methane at MV#2. In contrast, methane at MV#3 plots above the mixing line because of its high C₁/C₂. Below the SMTZ at MV#3, the δ¹³C-DIC values are high (ca. +20‰), and the isotope separation between δ¹³C-CH₄ and δ¹³C-DIC (expressed as the C isotopic separation factor ε_C ≈ δ¹³C_DIC − δ¹³C_CH4) was >70‰ (Figure 13C). Despite the relatively ¹³C-enriched CH₄ at MV#3, the large ε_C value indicates hydrogenotrophic methanogenesis (Whiticar, 1999). The addition of methane produced by hydrogenotrophic methanogenesis in a high-δ¹³C-DIC source sediment to thermogenic methane derived from greater depths was observed at a mud volcano in the Kumano forearc basin (Ijiri et al., 2018a). Such “secondary methanogenesis” may have occurred in the MV#3 mud reservoir (Milkov, 2011; Milkov and Etiope, 2018). As at MV#1–3, thermogenic methane should have been supplied to MV#14. However, because MV#14 is inactive, the present supply of thermogenic methane from deep sediments should be limited, and methanogenesis in shallower sediments likely overprinted the thermogenic signals at MV#14. Consequently, we assume that thermogenic methane was supplied to all studied mud volcanoes during their active periods, indicating that the supply of thermogenic hydrocarbon gases from a deep and hot sedimentary layer is characteristic of mud volcanism off Tanegashima Island. At IODP Site C0002 in the Kumano fore arc basin of the Nankai Trough, drilling to ∼2 kmbsf reached temperatures of ∼81°C–85°C, but the mixed thermogenic gas at that horizon was delivered from deeper and thus hotter sediments, indicating that no in situ thermogenic gas production occurs at that level (Hammerschmidt et al., 2014). Therefore, we infer that thermogenic gas production below the studied mud volcanoes does not occur below 80°C.