Upper-plate conduits linked to plate boundary that hosts slow earthquakes

Domain: earth_science
Published: 2023-09-20
Editorial outcome: minor_revision

Abstract
In shallow subduction zones, fluid behavior impacts various geodynamic processes capable of regulating slip behaviors and forming mud volcanoes. However, evidence of structures that control the fluid transfer within an overriding plate is limited and the physical properties at the source faults of slow earthquakes are poorly understood. Here we present high-resolution seismic velocity models and reflection images of the Hyuga-nada area, Japan, where the Kyushu-Palau ridge subducts. We image distinct kilometer-wide columns in the upper plate with reduced velocities that extend vertically from the seafloor down to 10–13 km depth. We interpret the low-velocity columns as damaged zones caused by seamount subduction and suggest that they serve as conduits, facilitating vertical fluid migration from the plate boundary. The lateral variation in upper-plate velocity and seismic reflectivity along the plate boundary correlates with the distribution of slow earthquakes, indicating that the upper-plate drainage system controls the complex pattern of seismic slip at subduction faults.

Introduction
Plate boundaries in subduction zones portray a wide spectrum of slip behavior, from day-to-year-long slow slip events to instantaneous faulting, governed by diverse physical factors. Pore fluid is one of the most influential parameters of seismic slips, and pore fluid behavior has been intensively studied in recent years1–3. High fluid pressure is a prerequisite condition for slow earthquakes to nucleate at subduction faults4. At shallow parts of subduction zones (depth < 15 km), the pore fluid pressure is primarily controlled by the compaction/dehydration reactions of sediments entering the subduction zone3. Once released from the source, the fluid migrates upward via permeable zones into the forearc wedge owing to its buoyancy5,6, forming mud volcanoes and cold seeps on the forearc basins and trench slopes7. Subducting seamounts (topographic highs at the plate boundary) impact the hydrological system in subduction zones by causing a heterogeneous distribution of permeability in the upper plate and forming a subduction channel filled with fluid-rich sediments8–10. However, the evidence on fluid distribution directly linked to subducting seamounts is limited. In particular, the along-strike extent of the impact of the subducting seamount is poorly constrained. In addition, our understanding of how the complex structure formed by subducting seamounts modulates the slip pattern at subduction faults remains unclear11.

The Hyuga-nada subduction zone in southwest Japan is an excellent area to study the relationship between subduction earthquakes and forearc hydrology. Adjacent to the Nankai Trough, where megathrust earthquakes of magnitudes greater than 8 occur frequently, the Hyuga-nada area historically hosted several interplate earthquakes of magnitudes ~712,13, the source areas of which are delineated by long-term slow slip events on the downdip side14 (Fig. 1a). On its updip side, slow earthquakes have been well-documented in recent seafloor observations15–18. For example, the activity of low-frequency tremors exhibiting rapid epicentral migration (with a speed of tens to hundreds of kilometers per day)15,16 implies that these events occur at the plate interface, driven by the spatiotemporal variation in the fluid pressure therein19. Another fluid-related phenomenon in the Hyuga-nada area is abundant mud volcanoes on the seafloor20–22. Although geochemical analyses of core samples suggested that these mud volcanoes are fed by the dehydration of clay minerals at shallow crustal depths21, the origin of the fluid is poorly understood. Moreover, the hydrological system responsible for the transfer of deep fluids is unknown.

To better understand the geodynamic system in the Hyuga-nada subduction zone, we performed a dense seismic refraction experiment using ocean bottom seismographs (OBSs) and a multichannel seismic (MCS) reflection survey in the summer of 2020. For the refraction study, we laid out 50 OBSs at 2-km intervals on the seismic line (HYU01) (Fig. 1b). The seismic line, aligned parallel to the regional trend of the Nankai Trough in the SW-NE direction, crossed the source area of the low-frequency tremor distributed between the shallow and deep seamounts (Fig. 1). In this study, we construct a P-wave velocity (Vp) model by applying the advanced full waveform inversion (FWI) technique to the OBS data (Methods section). Our results provide important structural constraints of the crust from the seafloor down to the subducting plate at a significantly higher spatial resolution than previous studies23–25. We compare the Vp model with the reflection image on the coincident line for structural interpretations. Based on the seismic results, we suggest the origin of deep fluids and their pathways up to the seafloor mud volcanoes. We also compare them with the spatial distribution of slow and regular earthquakes to discuss the structural controls on their occurrence. The structural characterization on the subduction megathrust faults using our study method can provide important information to evaluate the spatial extent of future megathrust earthquakes in the region.

Results
Upper-plate heterogeneity and plate boundary structure
The FWI Vp modeling along the HYU01 profile reveals a remarkably heterogeneous structure in the upper plate (Fig. 2a). The most important finding is the lateral variation in the Vp, highlighted by horizontally alternating high- and low-velocity anomalies (Fig. 2b), which is supported by the varying apparent velocities of first-arrival refraction phases in the OBS record sections (Supplementary Figs. 1–8). We detected several low-velocity columns having Vp values < 3.5 km/s that extended near-vertically from the seafloor to depths of over 10 km. These upper-plate velocities are significantly lower than those in other areas of the Nankai subduction zone located at similar distances from the trench axis (~50 km)25,26, suggesting that the subduction of the Kyushu-Palau Ridge increases the porosity of the forearc wedge over a broad region (about four times broader than the ridge itself). Interestingly, some of the low-velocity columns dip slightly toward the summit of the ridge located in the middle of the study area. Based on their fan-shaped geometry radiating from the seamount, we interpret that the low-velocity features represent an intensive fracture network; previous studies utilizing sandbox experiments explain that such networks consisting of subvertical faults can develop above subducting seamounts27. The dense distribution of fracture zones in the vicinity of the Kyushu-Palau Ridge may reflect the stress concentration in the region caused by seamount subduction, and their lower seismic velocities in the western part (compared with those in the eastern part; Fig. 2d–h) may indicate higher porosity and be attributed to the regional stress regime governed by the northwestward convergence of the incoming Philippine Sea Plate28.

In enlarged views of the seismic reflection image, the upper-plate fracture zones form conical zones with atypical reflection patterns characterized by inverted-V-shaped reflectors (Fig. 2c). The horizontal reflectors corresponding to the sedimentary layers are highly disturbed and portray upwardly domed geometry (Fig. 2i–m) on the top of the low-velocity fracture zones (Fig. 2d–h). These seismic reflection features are typical of mud volcanoes29 and suggest the occurrence of active upwelling flows inside the low-velocity columns. The seafloor bathymetry of the area shows several circular mounds (with steep slopes) close to the seismic profile, indicating another typical feature of mud volcanoes (Fig. 3). Additionally, the data suggest an extremely complex pattern of fault development, characterized by several linear fault scarps. Based on the consistency with experimental studies27, we suggest that these bathymetric features are associated with the ridge subduction. Importantly, the upper-plate low-velocity columns are situated where these linear scarps are developed (Fig. 3). These observational facts demonstrate that the upper plate accommodates well-developed fluid conduits that efficiently facilitate fluid flow to the seafloor and feed the mud volcanoes.

The seismic reflection image depicts reflective zones that extend laterally at depths of 10–13 km in the middle and eastern part of the seismic profile (40–100 km distance in Fig. 2c), with higher reflectivity in the flanks of the subducting ridge (40–45 km distance and 55–72 km distance in Fig. 2c) and away from the ridge (90–100 km distance in Fig. 2c). The reflective zones mostly trace the iso-Vp contour of 3.5–4.5 km/s (Fig. 2a) and also correspond to the areas having a positive velocity gradient in the depth direction, underlying those having a negative velocity gradient (Fig. 4). These layer structures probably represent low-velocity underthrust sediments atop the upper crust of the subducted Kyushu-Palau Ridge. In general, spatial variations in seismic reflectivity at the plate boundary is closely correlated with the physical properties and plate coupling conditions, which are highly influenced by in-situ fluid pressure30,31. At shallow depths where the temperature is below 150 °C, the pore fluid pressure is primarily controlled by the compaction and dehydration reactions of clay-rich sediments3. Notably, Vp values of 3.5–4.5 km/s imply that the reflective zones may be comprised of volcaniclastic sediments, which can hold a substantial amount of fluid before subduction and expel it when subducted10,32,33. Faint seismic reflections are observed in the western portion of the plate boundary. The disappearance of the reflection may suggest that fluids are starved there. Correspondingly, the westernmost low-velocity column is limited to the upper 8 km and unlikely to reach the plate boundary. The portion between the upper-plate low-velocity zones and the plate boundary increases the Vp to over 4.5 km/s and exhibits a sporadic reflection pattern with slightly higher reflectivity (0–40 km distance and 10–13 km depth in Fig. 2a and c).

Fluid origin and migration system for mud volcano formation
Submarine mud volcanoes are a well-known feature in forearc basins and provide a window to explore the subsurface hydrological system in subduction zones34. One of the key issues regarding mud volcanoes is the diversity of the fluid origins. Even though previous studies state that the fluids originate from the dehydration of clay-rich sediments within the accretionary prisms, which typically occurs at depths of several kilometer by rapid sedimentation and compaction34, recent studies report other potential sources, such as subducted sediments along deep thrust faults35 and forearc mantle wedges36. The seismic data recorded in our study offer new evidence for the previously derived conclusions that the fluid source for the mud volcanoes in the Hyuga-nada subduction zone is along the plate boundary at depths > 10 km.

At the shallow parts of subduction zones, fluid is likely to migrate via the decollement and branching faults, because they act as the most permeable pathways for fluid flow. A significant difference in the fluid feeding system in the Hyuga-nada area from this conventional concept is that the fluids may be ascending upward through the near-vertical conduits. We suggest that the fracture network in the Hyuga-nada subduction zone is primarily paved by the subducting seamount; the Vp model indicates that the network consists of discrete low-velocity columns that are densely distributed atop and around the subducting ridge (Fig. 2a and b). Our results also reveal some velocity reversals at different depths within the upper plate (Fig. 4). These results suggest that fluids from the plate boundary may be captured and accumulate at multiple depths within the upper plate and/or the upper plate itself also generates some fluids by dehydration of clay-rich sediments37. The latter possibility is also consistent with the observation that the upper-plate conduits are linked to not only the fluid-rich plate boundary (40–72 km and 80–100 km distances in Fig. 5a) but also the fluid-starved portions (30–40 km and 72–80 km distance in Fig. 5a).

Drainage system and its control on slow earthquakes
The variations in the Vp and reflectivity along the plate boundary agree well with the distribution of tremors and very low-frequency earthquakes (VLFEs) in the region. For example, while tremors and VLFEs actively occur around the subducting Kyushu-Palau Ridge, where the plate boundary is reflective and exhibits low seismic velocities (Vp of 3.5–4.5 km/s), their activities disappear in the western part of the profile (0–40 km distance in Fig. 5b), corresponding to the increase in the Vp value (> 4.5 km/s) and the waning seismic reflectivity around the plate boundary (0–40 km distance in Fig. 5a). Instead, the western part of the plate boundary is dominated by repeating regular earthquakes38 (Figs. 1 and 5b). This indicates a rapid transition in the slip pattern along the plate boundary that corresponds well with the variations in the physical properties around the plate boundary faults. As the Vp value measured in this study at the transitional area (~4.5 km/s) is consistent with that in previous studies on other subduction zones39,40, the value of 4.5 km/s is considered a common threshold to characterize the slip behavior along the megathrust faults.

Several episodes of tremors and VLFEs are observed within a localized area of an ~30-km range15–17 (“Isolated” tremor/VLFE episodes at 40–70 km distance in Fig. 5) and, similar to the western limit of the tremor zone, the eastern limits of these episodes coincide well with the nonreflective portion of the plate boundary (72–80 km distance in Fig. 2c). On the other hand, this nonreflective portion sometimes hosts the tremor activities (“Pervasive” tremor episodes in Fig. 5). Although factors controlling this variation in tremor occurrence are unclear, one possibility inferred from our structural observation is the role of the upper-plate conduit. The Vp model indicates that one of the low-velocity anomalies exits in the upper plate just above the non-reflective portion (at 72–80 km distance in Figs. 2b and 5a) and thus may suggest that the fluid pressure at the plate boundary is modulated by the partial leakage of fluids into the upper-plate conduits. In general, the upper-plate drainage system is known to play such a role4. Although this hypothesis is proposed for the tremor activity in deep subduction zones (i.e., depths > 30 km), it may also be true for shallow subduction zones (i.e., depths < 15 km).

Methods
Seismic data acquisition
The seismic data were acquired through two research cruises conducted in August 2020. During the KM20-05 cruise by the R/V Kaimei of Japan Agency for Marine-Earth Science and Technology (JAMSTEC), we deployed 50 OBSs along the HYU01 profile at 2-km intervals. Each OBS contained three‐component geophone and hydrophone sensors (with a natural frequency of 4.5 Hz; the sampling rate for recording was 200 Hz). The OBSs recorded the acoustic signals from 596 airgun shots triggered every 200 m along the 120-km-long seismic line. All the OBSs were recovered by the R/V Kairei during the KR20-10 cruise, providing the seismic records of the airgun signals. We also collected the MCS reflection data on the same line as the OBS refraction survey, using a 5.6-km-long 444-channel streamer cable onboard the Kaimei; the shooting interval for the reflection survey was 50 m. For both the surveys, the Kaimei’s tuned airgun arrays (with a total air volume of 10,600 cubic inch) were used as the source.

Vp model construction from OBS refraction data
After the recovery, the OBS data were preprocessed for further analyses. First, we located the OBS positions at the seafloor, using the same method applied in our previous study41. Then, we determined the first breaks in the refraction waves manually, using the vertical component of the geophone data. The picking error ranging from 30 ms to 60 ms were assigned to each pick, depending on the offset and the noise level of the data. A total of 27,902 first arrivals were used for a travel time tomography analysis42 to derive a preferable Vp model. The initial Vp model for the first-arrival travel time tomography (FAT) was created by referring to the seismic refraction studies close to the survey line considered in this study25 (Supplementary Fig. 9a). The final Vp model developed using FAT that explained the travel time data within the error range (root mean square [RMS] travel time residual of 56 ms) was obtained after 15 iterations of inversion (Supplementary Fig. 9b).

To improve the resolution of the Vp model, we further applied the FWI analysis to the OBS records. For the analysis, we used TOY2DAC, an acoustic 2D and frequency-domain FWI code developed by SEISCOPE26,43. For the modeling, the Vp structure of 100 × 23 km was gridded every 40 m. The result of FAT was used in the starting model. We applied the pre-processing to the vertical geophone data, such as muting before first arrival and time-damping with different time windows at each step of inversion. To stabilize the calculation, the inversion was initiated from a low frequency band of 2.5–3.5 Hz at the first step and then sequentially increased up to a higher frequency band of 2.5–7.5 Hz at the final step. We sampled the frequencies every 0.25 Hz within each frequency band. Other inversion parameters are summarized in Supplementary Table 1. The OBS records contained significant attenuations in the refraction phases, with shadow zones at some locations. To better reproduce these features in the OBS records, we introduced Vp-dependent attenuation factors (Q values) in the calculation. We used the Q values of 10,000, 50, 90, 350, 500, and 850 for the regions having Vp values < 1.55, 1.55–3.5, 3.5–4.5, 4.5–6.5, 6.5–7.5, and > 7.5 km/s, respectively44,45. The initial density model was converted from the Vp model, using the Gardner relation46. Both the Q and density values remained unchanged during inversion, and only the Vp values were updated. The overall waveform fits (Supplementary Fig. 1–8), reduction in the cost function (Supplementary Table 1), and consistency with the MCS reflection data were carefully examined to evaluate the reliability of the final model. The checkerboard resolution tests demonstrated that the FWI analyses was capable for resolving the velocity anomalies with 1.5 km (horizontally) x 1 km (vertically) down to the depth of ~15 km (Supplementary Fig. 10).

MCS reflection data processing
The primary processing of the MCS reflection data was carried out by DownUnder GeoSolutions (DUG). The data were processed following a standard prestack flow that included trace editing, tidal static correction, and low-cut filtering (of 1–3 Hz). The processing also included swell-noise attenuation, deghosting, source signature deconvolution, 2D surface-related multiple elimination, least-squares adaptive subtraction, and high-resolution radon demultiple for noise reduction. Velocity analyses and prestack Kirchhoff depth migration were performed to produce the final seismic reflection image (Fig. 2c).

Magnetic data processing
The magnetic anomaly data were downloaded from the Marine Trackline Geophysical Database of the National Centers for Environmental Information (https://www.ngdc.noaa.gov/mgg/gdas/). The data were gridded at 1 arc-min (bin size), and then were transformed to the vertical component of the magnetic field produced by the same source magnetized in the vertical direction [reduced-to-pole (RTP) anomaly]. The paleolatitude of the study area was assumed to be 26° N, while considering the northward convergence direction of the Philippine Sea Plate.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents a high-resolution two-dimensional seismic velocity model along the forearc wedge of the Hyuga-nada subduction zone obtained from full waveform inversion (FWI) of wide-angle ocean bottom seismic data. The model illuminates sub-vertical low velocity zones (LVZs) that the authors interpret as damaged and fractured zones in the upper plate created by subduction of a seamount that form a fluid drainage network that controls the hydrogeology system in this region. By comparing the seismic model with seafloor bathymetry, magnetic data, and non-volcanic tremor and very-low-frequency earthquakes, the authors make a convincing case that the LVZs represent a fluid drainage zones that impact the plate interface seismogenic properties.\nUnderstanding fluid migration and retention in subduction settings is topic of the greatest importance because of the major role that fluids play in subduction plate interface properties, with implication for seismic hazards. In this context, this manuscript will be an important contribution in this field and it would be an appropriate publication for Nature Communications.\nThe manuscript is succinct, well written, with nice and clear illustrations, and makes a convincing case for most of the interpretations. However there are some aspects that require improvement before acceptance, in particular more details about the FWI modeling approach, and discussion about some of the interpretations that are apparently not fully supported by the data.\nBelow is a list of all my comments, both major, minor, and editorial, in order of appearance in the manuscript. Numbers in my comments refer to numbered manuscript lines in the submitted PDF document.\nLines 45-46. Suggested rewording: “We image distinct several-kilometer-wide columns in the upper plate with reduced velocities that extend vertically from the seafloor up to 10-13 km depth.”.\n47. Change “radiated” to “radiate”, and “were” to “are”.\n51. Change “correlates to” to “correlate with”.\n66. Change “high” to “highs”.\n75. Delete “huge”.\n99. Change “the identical” to “a coincident”.\n114. Change “near-vertically to the depths” to “near-vertically from the seafloor to depths”.\n116. Change “that have the same distance from” to “located at similar distances from”.\n117-118. Suggested rewording: “…Kyushu-Palau Ridge increases the porosity of the forearc wedge over a broad region (about four times broader than the ridge itself).”.\n118. Change “incline” to “dip”.\n123. Delete “the” before fracture, and change “just above” to “in the vicinity”.\n125. Please define how “maturity” is used here. How does Vp or geometry of fractured zones inform about maturity?\n132. Change “which indicate a great correlation with” to “which correlate spatially with”.\n145. This sentence indicates that reflectivity is higher at ~40-60 km model distance. However reflectivity over the subducted seamount (~45-55 km) is actually lower than in other regions. Reflectivity is stronger away from the ridge at ~64-72 km and 90-100 km.\n149. Accreted sediments?\n150. Delete “the” preceding “spatial” and “seismic”. Change “variation” to “variations”.\n151. Change “correlated to” to “correlated with”.\n155, 156. Change “reflective zones are” to “reflective zones may be”, and “sediments that can” to “sediments, which can”.\n159-162. These sentences seem to imply that low plate reflectivity in the western section indicates fluid starvation, and it is reflected by LVZs being limited to the upper section. However the LVZ near 30 km model distance is underlaid by a narrow vertical LVZ extending to 12 km depth where plate boundary reflectivity is low. This structure seems to contradict the interpretation given in these sentences.\n172. Change “offer a new implication on” to “offer new evidence for”.\n181. Change “consisted of independent” to “consists of discrete”.\n206-210. The wording in these sentences is confusing. Please indicate which LV anomaly exits in the upper plate above the non-reflective interface (the anomaly shown in Fig. 2e, I presume?). Are the authors implying that efficient drainage of fluids through this anomaly reduces pore pressure and creates a barrier to tremor migration? If that is what is being suggested here, then it fails to explain why the anomaly at ~70-80 km, which is draining fluids from a non-reflective section of the plate interface (fluid starved, accordingly to Fig. 5a), does not act as a barrier to tremor in this section of the profile.\n258. Change “were applied to the traveltime” to “were used for a traveltime”.\n259. Reference 45 describes a traveltime mapping method (i.e., a wide-angle pre-stack depth migration), not a traveltime tomography method. Please give details or cite reference for the FAT analysis.\n269. Please indicate which data were used for the FWI: vertical geophone or hydrophone? Please also indicate it in the captions of Extended Data Figures 1-8.\n269-270. Please show in Extended Data Figures 1-8 the windows and mutes used so readers can see which portions of the record sections were used in the FWI.\n271. Are 2.5 and 7.5 Hz the only two frequencies used? Authors should show a graph of the data spectra.\n275-276. Based on the Vp-Q values listed here, the authors have used Q=90 for Vp=3.5-4.5 km/s. However in cited reference 47, Q=349 for Vp=3.6-4.6 km/s. Please explain the discrepancy, and comment on any potential implications for the FWI results.\n279. Change “inverted.” to “updated.”. Authors should also show a figure illustrating the reduce in cost function.\n293. Please include the pre-stack depth migration velocity model in Extended Data figure 9.\nFigure 5a. In the caption, indicate that background color is Vp as in Fig. 1. The slab Moho dashed line shown in this figure has not been shown in previous figures nor mentioned in the text. Cite reference for the source of this Moho interface.", "reviewer_2": "This manuscript connects the occurrence of slow earthquakes with the presence of fluids in the upper plate. Some structures observed in the seismic image are interpreted as fluid migration paths based on the correlation with low-velocity zones observed in the velocity model and with the location of mud volcanoes. The zones interpreted as fluid-rich broadly correlate to the spatial distribution of tremors and very low-frequency earthquakes.\nThe manuscript is well-written and explanatory. However, the influence of the fluid distribution in the seismogenic behaviour of subduction zones is already known and it is described in different publications (see a review in Saffer and Tobin, 2011; and Philip et al., 2023 for an example of a recent study in the Cascadia subduction zone). The changes in the reflectivity of the plate boundary reflector on seismic images were described in Ranero et al. (2008), and this reference should be included in the manuscript. Thus, the main novelty of these results may be the direct comparison with the recorded seismicity in the area. However, the comparison of the seismic images with the earthquake location is only done spatially. For example, in Figure 5c, earthquakes within a large area are included, but the seismic profile shows high variability in the structure of the fluid migration system in distances <15 km. This variability in the structure should be integrated into the discussion of the seismogenic behaviour. To strengthen the results, I suggest focusing on the earthquakes occurring close to the seismic profile, and locating these earthquakes along the profile, to compare the crustal structure with the seismicity of the area. Also, nothing is mentioned about the depth of nucleation of these earthquakes. Are they nucleated at the plate boundary or at upper crust depths? If they are located at upper crust depths, why is there a difference driven by the fluid content of the plate boundary?\nRegarding the fluid’s interpretation, there are fluids interpreted in the upper crust of the entire seismic profile, although there are only interpreted at the plate boundary towards the NE. The correlation of the low-velocity zones observed in the velocity model with mud volcanoes is unclear to me (Figure 3, see comments below). In my opinion, further discussion of the fluid origin and its implication for the seismogenic behaviour is needed.\nRegarding the figures, I suggest revising Figure 3 to make it easier to understand. As it is, the relation of the low-velocity column with the mud volcanoes is difficult to see. Are there any mud volcanoes at the exact location of the low-velocity columns?\nIn Figure 1 a, I can’t find the green shaded area.\nSuggested bibliography:\nPhilip, B. T., Solomon, E. A., Kelley, D. S., Tréhu, A. M., Whorley, T. L., Roland, E., Tominaga, M., & Collier, R. W. (2023). Fluid sources and overpressures within the central Cascadia Subduction Zone revealed by a warm, high-flux seafloor seep. Science Advances, 9(4). https://doi.org/10.1126/sciadv.add6688\nRanero, C. R., Grevemeyer, I., Sahling, H., Barckhausen, U., Hensen, C., Wallmann, K., Weinrebe, W., Vannucchi, P., Von Huene, R., & McIntosh, K. (2008). Hydrogeological system of erosional convergent margins and its influence on tectonics and interplate seismogenesis. Geochemistry, Geophysics, Geosystems, 9(3). https://doi.org/10.1029/2007GC001679\nSaffer, D. M., & Tobin, H. J. (2011). Hydrogeology and mechanics of subduction zone forearcs: Fluid flow and pore pressure. Annual Review of Earth and Planetary Sciences, 39, 157–186. https://doi.org/10.1146/annurev-earth-040610-133408"}, "author_response": "Point by point response to reviewer comments. Our responses are in red bold.\nReviewer #1\nThis manuscript presents a high-resolution two-dimensional seismic velocity model\nalong the forearc wedge of the Hyuga-nada subduction zone obtained from full\nwaveform inversion (FWI) of wide-angle ocean bottom seismic data. The model\nilluminates sub-vertical low velocity zones (LVZs) that the authors interpret as damaged\nand fractured zones in the upper plate created by subduction of a seamount that form a\nfluid drainage network that controls the hydrogeology system in this region. By\ncomparing the seismic model with seafloor bathymetry, magnetic data, and\nnon-volcanic tremor and very-low-frequency earthquakes, the authors make a\nconvincing case that the LVZs represent a fluid drainage zones that impact the plate\ninterface seismogenic properties.\nUnderstanding fluid migration and retention in subduction settings is topic of the\ngreatest importance because of the major role that fluids play in subduction plate\ninterface properties, with implication for seismic hazards. In this context, this\nmanuscript will be an important contribution in this field and it would be an appropriate\npublication for Nature Communications.\nThe manuscript is succinct, well written, with nice and clear illustrations, and makes a\nconvincing case for most of the interpretations. However there are some aspects that\nrequire improvement before acceptance, in particular more details about the FWI\nmodeling approach, and discussion about some of the interpretations that are apparently\nnot fully supported by the data.\nBelow is a list of all my comments, both major, minor, and editorial, in order of\nappearance in the manuscript. Numbers in my comments refer to numbered manuscript\nlines in the submitted PDF document.\nLines 45-46. Suggested rewording: “We image distinct several-kilometer-wide columns\nin the upper plate with reduced velocities that extend vertically from the seafloor up to\n10-13 km depth.”.\nWe modified as suggested (Lines 45-47 in the revised manuscript with track\nchanges).\n47. Change “radiated” to “radiate”, and “were” to “are”.\nWe deleted this part in the process of reducing the volume of abstract.\n51. Change “correlates to” to “correlate with”.\nWe replaced “to” with “with” (Line 53 in the revised manuscript with track\nchanges).\n66. Change “high” to “highs”.\nWe modified as suggested (Line 68 in the revised manuscript with track changes).\n75. Delete “huge”.\nWe modified as suggested (Line 77 in the revised manuscript with track changes).\n99. Change “the identical” to “a coincident”.\nWe replaced “identical” with “coincident” (Line 101 in the revised manuscript\nwith track changes).\n114. Change “near-vertically to the depths” to “near-vertically from the seafloor to\ndepths”.\nWe modified as suggested (Line 116 in the revised manuscript with track changes).\n116. Change “that have the same distance from” to “located at similar distances from”.\nWe modified as suggested (Line 118 in the revised manuscript with track changes).\n117-118. Suggested rewording: “…Kyushu-Palau Ridge increases the porosity of the\nforearc wedge over a broad region (about four times broader than the ridge itself).”.\nWe modified as suggested (Lines 120-121 in the revised manuscript with track\nchanges).\n118. Change “incline” to “dip”.\nWe modified as suggested (Line 122 in the revised manuscript with track changes).\n123. Delete “the” before fracture, and change “just above” to “in the vicinity”.\nWe modified as suggested (Line 127 in the revised manuscript with track changes).\n125. Please define how “maturity” is used here. How does Vp or geometry of fractured\nzones inform about maturity?\nWe agree that the term “maturity” was ambiguous here. We intended to suggest\nthat the western part has lower Vp values and thus higher porosity than eastern\npart. We modified in that way (Lines 129-131 in the revised manuscript with track\nchanges).\n132. Change “which indicate a great correlation with” to “which correlate spatially\nwith”.\nWe modified as suggested (Lines 137-138 in the revised manuscript with track\nchanges).\n145. This sentence indicates that reflectivity is higher at ~40-60 km model distance.\nHowever reflectivity over the subducted seamount (~45-55 km) is actually lower than in\nother regions. Reflectivity is stronger away from the ridge at ~64-72 km and 90-100 km.\nThank you for this comment. We revised the explanation accordingly (Lines\n150-152 in the revised manuscript with track changes).\n149. Accreted sediments?\nWe added “underthrust” before “sediments (Line 155 in the revised manuscript\nwith track changes).\n150. Delete “the” preceding “spatial” and “seismic”. Change “variation” to “variations”.\nWe modified as suggested (Lines 156-157 in the revised manuscript with track\nchanges).\n151. Change “correlated to” to “correlated with”.\nWe modified as suggested (Line 157 in the revised manuscript with track changes).\n155, 156. Change “reflective zones are” to “reflective zones may be”, and “sediments\nthat can” to “sediments, which can”.\nWe modified as suggested (Line 162 in the revised manuscript with track changes).\n159-162. These sentences seem to imply that low plate reflectivity in the western\nsection indicates fluid starvation, and it is reflected by LVZs being limited to the upper\nsection. However the LVZ near 30 km model distance is underlaid by a narrow vertical\nLVZ extending to 12 km depth where plate boundary reflectivity is low. This structure\nseems to contradict the interpretation given in these sentences.\n“the westernmost low-velocity column” (Line 167 in the revised manuscript with\ntrack changes) indicates not the LVZ near 30 km model distance but one at 15-20\nkm model distance. A body with Vp over 4.5 km/s obviously underlies the LVZ at\n15-20 km model distance at depths greater than 8 km and thus the LVZ is unlikely\nto reach the plate boundary as described in the text (Lines 167-168 in the revised\nmanuscript with track changes). It is unclear whether the LVZ near 30 km model\ndistance reach the plate boundary since the same Vp body exists at the bottom of\nthe LVZ (See our interpretations in Figure 5a).\n172. Change “offer a new implication on” to “offer new evidence for”.\nWe modified as suggested (Lines 178-179 in the revised manuscript with track\nchanges).\n181. Change “consisted of independent” to “consists of discrete”.\nWe modified as suggested (Line 187 in the revised manuscript with track changes).\n206-210. The wording in these sentences is confusing. Please indicate which LV\nanomaly exits in the upper plate above the non-reflective interface (the anomaly shown\nin Fig. 2e, I presume?). Are the authors implying that efficient drainage of fluids\nthrough this anomaly reduces pore pressure and creates a barrier to tremor migration? If\nthat is what is being suggested here, then it fails to explain why the anomaly at ~70-80\nkm, which is draining fluids from a non-reflective section of the plate interface (fluid\nstarved, accordingly to Fig. 5a), does not act as a barrier to tremor in this section of the\nprofile.\nTo avoid this confusion, we revised the text to clearly indicate that the LV anomaly\nmentioned here is the one at 70-80 km distance (not the western one in Fig. 2e;\nLines 215-216 in the revised manuscript with track changes). As mentioned in the\ntext (Lines 209-213 in the revised manuscript with track changes) and shown in Fig.\n5, 70-80 km distance in the profile corresponds to the eastern end of some of the\ntremor episodes and thus seems to act as a barrier.\n258. Change “were applied to the traveltime” to “were used for a traveltime”.\nWe modified as suggested (Line 267 in the revised manuscript with track changes).\n259. Reference 45 describes a traveltime mapping method (i.e., a wide-angle pre-stack\ndepth migration), not a traveltime tomography method. Please give details or cite\nreference for the FAT analysis.\nWe replaced reference 45 with Fujie et al. (2013) (new reference 46; Lines 504-509\nin the revised manuscript with track changes).\n269. Please indicate which data were used for the FWI: vertical geophone or\nhydrophone? Please also indicate it in the captions of Extended Data Figures 1-8.\nWe used vertical geophone data. We added the description to Line 279 and\ncaptions of Supplementary Figures 1-8.\n269-270. Please show in Extended Data Figures 1-8 the windows and mutes used so\nreaders can see which portions of the record sections were used in the FWI.\nWe added an additional figure showing the offset range and the time windows used\nin the FWI to the Supplementary Information File (Supplementary Figure 11).\n271. Are 2.5 and 7.5 Hz the only two frequencies used? Authors should show a graph of\nthe data spectra.\nWe modified this statement as we used every 0.25 Hz between 2.5 and 7.5 Hz\n(Lines 281-285 in the revised manuscript with track changes). We also added a\ntable summarizing inversion parameters used at each inversion step\n(Supplementary Table 1).\n275-276. Based on the Vp-Q values listed here, the authors have used Q=90 for\nVp=3.5-4.5 km/s. However in cited reference 47, Q=349 for Vp=3.6-4.6 km/s. Please\nexplain the discrepancy, and comment on any potential implications for the FWI results.\nWe used the Q value of 90 for accretionary prisms from reference 48 (Table 2 in\nFurumura et al., 2008) because we confirmed that this low Q value can reproduce\nsynthetic waveforms that agree well with observed records. We also tested with a\nhigh Q value of 349 and found the inversion tend to be unstable very easily, so we\nderived only the Q value of 50 for the shallowest part of sediments from reference\n47 (Hino et al., 2015). Hino et al. (2015) discuss that Q value of 349 from the\nKumano Basin is much higher than other forearc areas, such as Costa Rica where\nQp is estimated to be 50 to 150 (Zhu et al., 2010). Hino et al. (2015) raise a\npotential cause of the difference that the Kumano Basin is rigid (relatively stable)\naccretionary prism while the Costa Rica is an erosive margin associated with\nseamount subduction and thus host extensive fracturing in the forearc. We think\nthe tectonic setting in the Hyuga-nada area is more similar as the case of Costa\nRica.\n279. Change “inverted.” to “updated.”. Authors should also show a figure illustrating\nthe reduce in cost function.\nWe changed “inverted” to “updated” as suggested (Line 292 in the revised\nmanuscript with track changes). For the reduction in cost function, it was difficult\nto illustrate it in a single figure or graph because we use different frequency values,\noffset range, and time windows at each inversion step. Instead of that, we describe\nit by the ratio of initial and final models at each inversion step together with other\ninversion parameters in Supplementary Table 1.\n293. Please include the pre-stack depth migration velocity model in Extended Data\nfigure 9.\nWe added the pre-stack depth migration velocity model to Supplementary Figure\n9.\nFigure 5a. In the caption, indicate that background color is Vp as in Fig. 1. The slab\nMoho dashed line shown in this figure has not been shown in previous figures nor\nmentioned in the text. Cite reference for the source of this Moho interface.\nWe added a statement to the caption of Figure 5a that indicate the background\ncolor is Vp in Figure 2a. We used iso-Vp contour of 7.5 km/s as a proxy of the slab\nMoho. We added this statement as well (Lines 578-579 in the revised manuscript\nwith track changes).\nThis manuscript connects the occurrence of slow earthquakes with the presence of fluids\nin the upper plate. Some structures observed in the seismic image are interpreted as\nfluid migration paths based on the correlation with low-velocity zones observed in the\nvelocity model and with the location of mud volcanoes. The zones interpreted as\nfluid-rich broadly correlate to the spatial distribution of tremors and very low-frequency\nearthquakes.\nThe manuscript is well-written and explanatory. However, the influence of the fluid\ndistribution in the seismogenic behaviour of subduction zones is already known and it is\ndescribed in different publications (see a review in Saffer and Tobin, 2011; and Philip et\nal., 2023 for an example of a recent study in the Cascadia subduction zone). The\nchanges in the reflectivity of the plate boundary reflector on seismic images were\ndescribed in Ranero et al. (2008), and this reference should be included in the\nmanuscript.\nOne of the novelties of this study is a discovery of vertical fluid conduits in the\nupper plate. This structure is obviously different from typical accretionary wedge\nsystem and, to our knowledge, has never been observed and discussed in any\npublications. In addition, we relate them to ridge subduction based on the\ncomparison with seafloor bathymetry and sandbox experiments. We understand\nthat many previous studies (including Ranero et al., 2008) have examined the\nrelation between the seismic reflectivity and seismogenesis at the plate interface in\nsubduction zones worldwide, but never been done in the Hyuga-nada area. Saffer\nand Tobin (2011) has been cited already (ref 3).\nThus, the main novelty of these results may be the direct comparison with the recorded\nseismicity in the area. However, the comparison of the seismic images with the\nearthquake location is only done spatially. For example, in Figure 5c, earthquakes\nwithin a large area are included, but the seismic profile shows high variability in the\nstructure of the fluid migration system in distances <15 km. This variability in the\nstructure should be integrated into the discussion of the seismogenic behaviour. To\nstrengthen the results, I suggest focusing on the earthquakes occurring close to the\nseismic profile, and locating these earthquakes along the profile, to compare the crustal\nstructure with the seismicity of the area. Also, nothing is mentioned about the depth of\nnucleation of these earthquakes. Are they nucleated at the plate boundary or at upper\ncrust depths? If they are located at upper crust depths, why is there a difference driven\nby the fluid content of the plate boundary?\nWe have already incorporated most of the suggestions here in the original\nmanuscript. To focus on the tremor and VLFE that occurred close to our seismic\nprofile, we only used the events within 5-km distance from the seismic profile (See\nthe caption of Figure 5b). Unfortunately, as the depth of these events are poorly\nconstrained as explained by the cited references (Yamashita et al., 2015, 2021;\nTonegawa et al., 2020), we cannot distinguish whether they occur at the plate\nboundary, in the upper plate or both. We do not suggest anywhere that the events\nare located in the upper plate. All we can do is to discuss their spatial correlation in\nthe horizontal directions as in the text.\nRegarding the fluid’s interpretation, there are fluids interpreted in the upper crust of the\nentire seismic profile, although there are only interpreted at the plate boundary towards\nthe NE. The correlation of the low-velocity zones observed in the velocity model with\nmud volcanoes is unclear to me (Figure 3, see comments below). In my opinion, further\ndiscussion of the fluid origin and its implication for the seismogenic behaviour is\nneeded.\nIn the southwestern part of the profile, the plate boundary is deeper than 15 km,\nwhich may be too deep for underthrust sediments to generate fluids via\ndehydration. To explain the slight difference in fluid distribution between the\nupper plate and the plate boundary, we added another possibility on the fluid\norigins, that is dehydration of clay-rich sediments within the upper plate, based on\na new reference (Mitsutome et al. 2023, Scientific Reports) (Lines 190-192 in the\nrevised manuscript with track changes).\nRegarding the figures, I suggest revising Figure 3 to make it easier to understand. As it\nis, the relation of the low-velocity column with the mud volcanoes is difficult to see.\nAre there any mud volcanoes at the exact location of the low-velocity columns?\nThe relationship between the upper-plate low-velocity columns and the mud\nvolcanoes on the seafloor is not straightforward because many of the low-velocity\ncolumns are still buried in the sedimentary layers (Fig. 2d-2i, 2e-2j, 2f-2k, 2g-2l).\nThere is only one example forming a clear mound structure on the seafloor (Fig.\n2h-2m). Figure 3 also aims to suggest that the study area hosts not only circular\nmounds indicating mud volcanoes but also a number of linear scarps associated\nwith the ridge subduction, which is consistent with sandbox experiments\n(Dominguez et al. 1998).\nIn Figure 1 a, I can’t find the green shaded area.\nTwo small green shaded areas are located just east of Kyushu Island.\nSuggested bibliography:\nPhilip, B. T., Solomon, E. A., Kelley, D. S., Tréhu, A. M., Whorley, T. L., Roland, E.,\nTominaga, M., & Collier, R. W. (2023). Fluid sources and overpressures within the\ncentral Cascadia Subduction Zone revealed by a warm, high-flux seafloor seep. Science\nAdvances, 9(4). https://doi.org/10.1126/sciadv.add6688\nRanero, C. R., Grevemeyer, I., Sahling, H., Barckhausen, U., Hensen, C., Wallmann, K.,\nWeinrebe, W., Vannucchi, P., Von Huene, R., & McIntosh, K. (2008). Hydrogeological\nsystem of erosional convergent margins and its influence on tectonics and interplate\nseismogenesis. Geochemistry, Geophysics, Geosystems, 9(3).\nhttps://doi.org/10.1029/2007GC001679\nSaffer, D. M., & Tobin, H. J. (2011). Hydrogeology and mechanics of subduction zone\nforearcs: Fluid flow and pore pressure. Annual Review of Earth and Planetary Sciences,\n39, 157–186. https://doi.org/10.1146/annurev-earth-040610-133408"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Dear Editor,\nthe authors of manuscript NCOMMS-23-09256A have adequately addressed all the concerns and comments I had in my first review of this manuscript. I have no further comments on the revised manuscript. My recommendation is to publish the manuscript in its present form. Thank you for the opportunity to review this manuscript.", "reviewer_2": "I appreciate the effort that the authors have made to revise their manuscript. My comments have been tackled in the response letter. However, I still have some concerns mainly regarding the discussion of the interpretation and the seismogenic behavior. I understand that the length of the manuscript is limited, but from my point of view, some of the discussions are too short to fully convince the reader about the proposed interpretation. Below, there is a list of my detailed comments:\n1) Given the high resolution of the data presented, I would expect a more detailed interpretation of the fluid migration paths and the correlation of their expression among the velocity model and the MCS image, as well as their relation to the bathymetric features. I do not necessarily disagree with the proposed interpretation in Figure 5a, but I think that further discussion is needed to fully justify the differences along the profile. For example, the sentence “ The seismic data recorded in our study offer new evidence for the previously derived conclusions that the fluid source for the mud volcanoes in the Hyuga-nada subduction zone is along the plate boundary at depths >10 km” (lines 178-180 manuscript with track changes) is not in full agreement with the Figure 5a interpretation, as some of the fluid migration paths are interpreted as coming from the fluid-starved plate boundary (70-80 km) and also, the fluids interpreted between 15-35 km are associated with a fluid-starved plate boundary area and are interpreted as intra-crustal generated fluids. I would appreciate a longer discussion about the reasons behind this interpretation and the possibilities of the fluid origin. I think this is needed to justify your interpretation and it will strengthen your conclusions, as one of the main objectives of this manuscript is to improve the knowledge about fluid origin and migration paths.\n2) I think that the high-Vp body in Figure 50 (0-30 km) is interpreted solely from the velocity model. Do you see any evidence of this body in the MCS? Is it somehow related to the fluid generation in this area?\n3) For me, Figure 3 is not helping. As a suggestion, it would be nice to show the location of the mound that reaches the seafloor (figure 2h-m) on the map to show the correlation between both features. Also, I do not think that from figure 3 the straightforward conclusion is that the linear features are scarps due to the seamount subduction. I suggest revising this figure to clarify it.\n4) Regarding the seismicity distribution and as a follow-up to my previous model about the earth-quake depth: A discussion on how the fluid distribution affects the seismicity is made. However, the distribution of fluid migration paths along the profile is very similar (see Figure 5a), and the main change is whether the plate boundary is fluid-rich or fluid-starved. However, the content of fluids at the plate boundary does not seem to be related to the presence or absence of fluid flows in the upper crust. Thus, I do not fully understand the discussion about the changes in the seismogenic behaviour in relation to the partial leakage of fluids into the upper plate. Also, the presence of the high-vp velocity body (0-35 km) may have an effect on the seismicity distribution, and this possibility is not mentioned/discussed in the manuscript. In my opinion, a few more sentences about these issues will be helpful to further justify and explain your conclusions."}, "author_response": "Point by point response to reviewer comments. Our responses are in red bold.\nReviewer #1 (Remarks to the Author):\nDear Editor,\nthe authors of manuscript NCOMMS-23-09256A have adequately addressed all the\nconcerns and comments I had in my first review of this manuscript. I have no further\ncomments on the revised manuscript. My recommendation is to publish the manuscript\nin its present form. Thank you for the opportunity to review this manuscript.\nReviewer #2 (Remarks to the Author):\nI appreciate the effort that the authors have made to revise their manuscript. My\ncomments have been tackled in the response letter. However, I still have some concerns\nmainly regarding the discussion of the interpretation and the seismogenic behavior. I\nunderstand that the length of the manuscript is limited, but from my point of view, some\nof the discussions are too short to fully convince the reader about the proposed\ninterpretation. Below, there is a list of my detailed comments:\n1) Given the high resolution of the data presented, I would expect a more detailed\ninterpretation of the fluid migration paths and the correlation of their expression among\nthe velocity model and the MCS image, as well as their relation to the bathymetric\nfeatures. I do not necessarily disagree with the proposed interpretation in Figure 5a, but\nI think that further discussion is needed to fully justify the differences along the profile.\nFor example, the sentence “ The seismic data recorded in our study offer new evidence\nfor the previously derived conclusions that the fluid source for the mud volcanoes in the\nHyuga-nada subduction zone is along the plate boundary at depths >10 km” (lines\n178-180 manuscript with track changes) is not in full agreement with the Figure 5a\ninterpretation, as some of the fluid migration paths are interpreted as coming from the\nfluid-starved plate boundary (70-80 km) and also, the fluids interpreted between 15-35\nkm are associated with a fluid-starved plate boundary area and are interpreted as\nintra-crustal generated fluids. I would appreciate a longer discussion about the reasons\nbehind this interpretation and the possibilities of the fluid origin. I think this is needed\nto justify your interpretation and it will strengthen your conclusions, as one of the main\nobjectives of this manuscript is to improve the knowledge about fluid origin and\nmigration paths.\nThank you for this suggestion. To respond to this comment, we added a sentence\ndescribing the relationship of the upper-plate conduits with the plate boundary\n(Lines 189-192 in the revised manuscript with track changes). For its relation to\nthe bathymetry, please see our response to the third comment below and the\nupdates of Figure 3.\n2) I think that the high-Vp body in Figure 50 (0-30 km) is interpreted solely from the\nvelocity model. Do you see any evidence of this body in the MCS? Is it somehow\nrelated to the fluid generation in this area?\nThank you for pointing this out. The high-Vp body shows relatively high\nreflectivity in the MCS image, but we are not sure if the change in reflectivity is\nrelated to fluid generation as it shows higher Vp. We only describe about the\nreflectivity in the text (Lines 157-165 in the revised manuscript with track\nchanges).\n3) For me, Figure 3 is not helping. As a suggestion, it would be nice to show the\nlocation of the mound that reaches the seafloor (figure 2h-m) on the map to show the\ncorrelation between both features. Also, I do not think that from figure 3 the\nstraightforward conclusion is that the linear features are scarps due to the seamount\nsubduction. I suggest revising this figure to clarify it.\nWe agree that the bathymetric features are so complicated that it is not easy to\nunderstand them only from the present form of Fig. 3. However, we believe that\nthe linear features are primarily caused by the ridge subduction because their\nspatial pattern is similar to what is predicted by sandbox experiments that examine\nseamount subduction (ref 27). For this, we added some explanations to the text\n(Lines 135-136 in the revised manuscript with track changes). For better\nunderstanding of the spatial relationship between the subducting seamounts,\nupper-plate low-velocity zones and seafloor mounts and scarps, we also prepared a\nshaded relief map with iso-depth contours in Fig 3a. We added the locations of Fig.\n2d-2h to the maps as suggested.\n4) Regarding the seismicity distribution and as a follow-up to my previous model about\nthe earth-quake depth: A discussion on how the fluid distribution affects the seismicity\nis made. However, the distribution of fluid migration paths along the profile is very\nsimilar (see Figure 5a), and the main change is whether the plate boundary is fluid-rich\nor fluid-starved. However, the content of fluids at the plate boundary does not seem to\nbe related to the presence or absence of fluid flows in the upper crust. Thus, I do not\nfully understand the discussion about the changes in the seismogenic behaviour in\nrelation to the partial leakage of fluids into the upper plate. Also, the presence of the\nhigh-vp velocity body (0-35 km) may have an effect on the seismicity distribution, and\nthis possibility is not mentioned/discussed in the manuscript. In my opinion, a few more\nsentences about these issues will be helpful to further justify and explain your\nconclusions.\nThank you also for this comment. We agree that the tremor distribution seems\nrelated only to the reflectivity (fluid contents) at the plate boundary and is not\nalways relevant to the upper-plate low-velocity zones. To clarify this complexity\nand avoid potential confusions, we modified the statements (Lines 213-224 in the\nrevised manuscript with track changes). Although the discussion on the correlation\nof high-Vp body in the western part (0-35 km distance) with the seismicity\ndistribution has been made in the original manuscript, we modified the text for\nclarity (Lines 199-203 in the revised manuscript with track changes)."}]

Decision letter