Shallow Low‐Velocity Layer in the Hyuga‐Nada Accretionary Prism and Its Hydrological Implications: Insights From a Passive Seismic Array

Shear wave velocity (Vs) estimations of accretionary prisms can pose unique constraints to the physical properties of rocks, which are hard to obtain from compressional velocities (Vp) alone. Thus, it would help better understand the fluid processes of the accretion system. This study investigates the Vs structure of the Hyuga‐nada accretionary prism using an array of ocean‐bottom seismometers (OBSs) with a 2 km radius. Teleseismic Green's functions and a surface wave dispersion curve are inverted to one‐dimensional Vs structures using transdimensional inversion. The results indicate the presence of a low‐velocity zone 3–4 km below the seafloor. The reduced Vs corresponds to a reduced Vp feature obtained in a previous seismic refraction survey, and the high Vp/Vs ratio suggests the presence of high pore fluid pressure. In addition, the resolved lithological boundary exhibits a sharp offset that extends laterally across the OBS array. We attribute this offset to a blind fault below while acknowledging other possibilities, such as due to mud diapirism or intense fracturing. The predicted fault is located at the Kyushu–Palau Ridge flank, oriented roughly parallel to the ridge axis, and thus likely caused by ridge subduction. This fault may act as a fluid conduit, contributing to the formation of a fluid reservoir beneath the compacted sediment layers.

. Tectonic setting of the study area and array configuration. (a) The orange star denotes the location of an array of ocean-bottom seismometers (OBSs). The red line represents the cross-section shown in (c) and (d). Yellow dots represent the epicenters of the tectonic tremors (Yamashita et al., 2015(Yamashita et al., , 2021. The pink line denotes the subducted Kyushu-Palau Ridge (Yamamoto et al., 2013). (b) Array configuration. The gray contour indicates the water depth, with an interval of 10 m. (c) P-wave velocity model obtained from a refraction survey (Nakanishi et al., 2018). The yellow inverse triangles represent the locations of OBSs. The subducting plate interface is denoted by the black line, which is defined by the velocity gradient profile of (d). (d) The same as (c), but vertical velocity gradients are shown. PHP: Philippine Sea Plate.
3 of 15 subduction accompanies the right-lateral motion (Itoh et al., 1998;Yamazaki & Okamura, 1989). Tectonic tremors and very-low-frequency earthquakes, both are members of slow earthquakes, intermittently occur near the KPR with an interval of 1-3 years (Baba et al., 2020;Tonegawa et al., 2020;Yamashita et al., 2015Yamashita et al., , 2021. As suggested for other regions worldwide, these slow earthquake activities may reflect a fluid-rich environment near the plate interface (Saffer & Wallace, 2015). However, little is known about the fluid processes (e.g., fluid sources, pathways, reservoirs) in this region.
High-resolution structures of the accretionary prism in this region were obtained in previous active-source seismic surveys (Nakanishi et al., 2018;Nishizawa et al., 2009;Park et al., 2009). Figure 1c shows a P-wave velocity (Vp) model based on a refraction survey (Nakanishi et al., 2018). Overall, the accretionary prism shows Vp of 2-4 km/s, and the subducting Philippine Sea Plate has a higher velocity of >6 km/s beneath the prism. Interestingly, velocity inversion with depth is noticeable at ∼2 km beneath the seafloor (Figure 1d). Nishizawa et al. (2009) reported a similar low-velocity zone (LVZ) beneath another independent seismic profile in Hyuga-nada. These LVZs may indicate fluid-rich conditions, although previous studies have not provided a detailed interpretation. The challenges are the modest sensitivity of the refraction surveys to thin LVZs with a sharp velocity contrast and the interpretation of physical properties based on Vp alone.
This study investigates the shear wave velocity (Vs) structure by utilizing a dense passive seismic array of ocean-bottom seismometers (OBSs) deployed in the Hyuga-nada region. Traditionally, active-source seismic surveys play a central role in constraining Vs structures within shallow marine sediments (e.g., Tsuji et al., 2011). However, in contrast to Vp, investigating Vs via active seismic sources is challenging because of the inefficient excitation of shear waves beneath the seafloor. In recent years, various elements of passive seismic records have been increasingly used to overcome this problem, including ambient surface wave noise (Mosher et al., 2021;Tonegawa et al., 2017;Yamaya et al., 2021;Zhang et al., 2020), teleseismic body waves (Agius et al., 2018;Akuhara et al., 2020), and a combination of them (Doran & Laske, 2019). This study attempts to solve Vs structures through the transdimensional inversion of teleseismic body waves and a surface wave dispersion curve (DC). Based on the results, we discuss the hydrological features in Hyuga-nada, which can be linked to the subducted KPR.

Passive Seismic Array
This study uses a passive seismic array of 10 OBSs installed in the Hyuga-nada region. The OBSs continuously recorded seismic waveforms from 30 March 2018 to 30 September 2018 ( Figure 1). Five OBSs (HDA01-05) were evenly installed within a radius of 1 km, whereas the other five OBSs (HDA06-10) were placed within 2 km, around the same center. Each OBS contains short-period three-component sensors (LE-3Dlite, Lennartz Electronic GmbH, Germany) and a gimbal to maintain the sensor's horizontality. The seismometer positions were constrained by acoustic positioning from a research vessel. The sensor orientations were determined from the particle motion of teleseismic Rayleigh waves (Sawaki et al., 2022;Stachnik et al., 2012; see text S1, Figure S1, and Table S1 in Supporting Information S1).
The array aimed to explore the potential of passive source methods for imaging shallow sediment structures. Another broadband OBS was deployed at the center of the array circle, but we failed to recover it. The array was placed on the refraction seismic survey line such that the tomography model could be used as a reference (Nakanishi et al., 2018; Figure 1). According to this refraction survey, the interface of the Philippine Sea Plate subducts to ∼10-11 km depth beneath the array. The seafloor topography is relatively gentle, with a slight slope to the northeast, resulting in a height difference of only ∼120 m over the 4 km diameter (Figure 1b). Therefore, its effects on surface and body wave propagation are negligible.

Method
This section elaborates on procedures we adopted to estimate Vs structures beneath the OBSs. The DC measurements from ambient noise records are described in Section 3.1. In Section 3.2, we describe the procedure we used to retrieve the Green's function (GF) from teleseismic P-waves. Subsequently, the acquired DC and GFs were inverted to one-dimensional (1D) Vs structures using a transdimensional, stochastic inversion scheme, as discussed in Section 3.3.

Rayleigh Wave Dispersion Curve
We retrieved Rayleigh waves propagating across the array from half-year-long records of ambient seismic noise. For this purpose, we employ a series of signal processing steps: cutting records into one-hour-long segments, detrending time series, downsampling data from 200 to 10 samples per second, deconvolution with instrumental responses, spectral whitening, and one-bit normalization in the time domain (Bensen et al., 2007). Cross-correlation functions (CCFs) of vertical components are then calculated between each station pair and stacked over the entire observation period. Figure 2 shows vertical-component CCFs obtained from all station pairs. The fundamental Rayleigh mode dominates CCFs at 0.2-0.4 Hz with an apparent velocity of ∼0.5 km/s. Based on the assumption of a laterally homogeneous structure beneath the array, the aggregation of CCFs in Figure 2 can be considered a virtual short gather recorded by a linear array. We estimate an averaged DC across the OBS array by applying the frequency-wavenumber (FK) analysis to this virtual shot gather. This treatment can significantly extend the high-frequency (or short-wavelength) limit of phase velocity measurements without suffering from spatial aliasing effects (Gouédard et al., 2008). This feature benefits this study because acquiring higher-frequency phase velocities is essential to constrain shallow structures of marine sediments.
The FK domain spectrum obtained from these virtual records shows the DC of the fundamental Rayleigh wave, which is traceable from 0.15 Hz (near the resolution limit) to 0.5 Hz (Figure 3). The DC shows minor deflection at ∼0.4 Hz, which we consider an artifact from the specific array configuration. In the higher frequency range between 0.5 and 1.0 Hz, the spectrum exhibits a complex pattern, and it is hard to distinguish the actual signal from artificial sidelobes. A relatively continuous feature can be observed at a frequency of >1 Hz, corresponding to the higher-mode Rayleigh wave, but the mode identification is nontrivial because of the ambiguity in the range of 0.5-1.0 Hz.

Teleseismic Green's Functions
We extract P waveforms of teleseismic events with M > 5.5 and an epicentral distance of 30-90°. Each extracted record is decimated to 20 samples per second, and two horizontal components were rotated to radial and transverse directions. We only retain data with a signal-to-noise ratio (SNR) above 3.0 on the vertical component. In this study, the SNR is defined as the root-mean-square amplitude ratio of 30 s time windows before and after P arrival. The GFs of teleseismic P-waves are retrieved from these time windows with the blind deconvolution technique (Akuhara et al., 2019). In contrast to conventional receiver function methods that only solve radial-component GFs, both radial-and vertical-component GFs can be estimated with this method. The retrieval of vertical-component GFs is crucial for ocean-bottom settings because intense water multiples dominate the vertical-component records. We use 60 s time windows for the deconvolution and apply a Gaussian low-pass filter to the results. The Gaussian parameter (i.e., standard deviation) is set to 8, corresponding to a 10% gain at ∼4 Hz.
The radial-component GFs are mostly coherent across the array, especially for the first 4 s (Figure 4a; see also Figures S2-S11 in Supporting Information S1 for wiggle plots against event back azimuths). A negative peak is predominant at ∼2.0-2.5 s after the direct P arrival. This coherency quantitatively justifies the 1D structure assumption we made for the FK analysis. At zero lag time, a peak corresponding to the direct P arrival is not evident, indicating the nearly vertical incidence of the P phase due to the low Vp   Figure 4b). The first reverberation with a positive polarity is evident at 3.1 s, and the second one can be observed at 6.2 s and has a reversed polarity. Although we did not use these vertical-component GFs for the inversion analysis, the good recovery of water reverberations to some degree validates the radial component estimations.

Transdimensional Bayesian Inversion
We use a transdimensional Bayesian interface and the reversible-jump Markov chain Monte Carlo (RJMCMC) algorithm (Green, 1995) for the inversion of the dispersion and GF data to an isotropic 1D Vs model beneath each OBS. The RJMCMC performs probabilistic sampling of model parameters, allowing the dimension of the model parameter space to be unknown. In our case, the algorithm automatically selects the number of layers in the 1D subsurface structure model. The transdimensional Bayesian inversion aims to estimate the posterior probability of the model parameter, m k , with the given data, d, that is, ( | ) , where k is a parameter determining the model-space dimension. Based on the Bayes' theorem, the posterior probability is proportional to the product of the prior probability, ( ) , and the likelihood, ( | ) :

Model Parameters
We assume that the subsurface structure consists of k layers. Each layer has constant seismic P-and S-wave velocities and density; the structure's lateral heterogeneity, anisotropy, and dissipation are ignored. We defined a model vector = ( 1, ⋯, −1, 1, ⋯, −1, DC, GF) , where δβ i is the S-wave velocity perturbation relative to a reference model and z i is the bottom depth of the ith layer. The other two parameters, σ DC and σ GF , represent the standard deviations of data noise, which are also solved within the hierarchical Bayesian model (Bodin et al., 2012). Based on a given set of model parameters, first, a Vs value of each layer is extracted from the reference model. The perturbation δβ i is then added to the extracted value. Similarly, Vp is obtained from the reference model, but without perturbation. The density is calculated from the Vp using an empirical relationship (Brocher, 2005). We fix the properties of the bottom half-space (i.e., k th layer) to stabilize the forward calculation of dispersion curves: Vs is set to 4.0 km/s and Vp and the density are scaled to Vs using the empirical law of Brocher (2005). For the seawater layer, we assume an acoustic velocity of 1.5 km/s and thickness of 2.388 km, which is the average station depth. The reference model was constructed from the two-dimensional (2D) P-wave velocity model of Nakanishi et al. (2018), as shown in Figure 1c, with the empirical scaling law that converts Vp into Vs (Brocher, 2005). Since the lateral velocity variation is small across the array, we construct a single reference model and apply it to all stations.

Likelihood
We calculate the likelihood ( | ) based on the assumption of Gaussian noise distribution: where C DC and C GF are the covariance matrices of the DC and GF data noise, respectively, and g DC and g GF are the synthetic DC and GF, respectively. The data vector, d, consists of DC and GF data vectors, denoted as d DC and d GF , respectively, with a length of N DC and N GF , respectively. We assume the temporal correlation of noise for GFs, which originates from the Gaussian low-pass filter, and a constant noise level across the entire time series. The corresponding covariance matrix can be expressed by CGF = 2 GF ( − ) 2 , where r is pre-determined from the Gaussian filter width (Bodin et al., 2012) and σ GF is a standard deviation of the data noise. We ignore off-diagonal components of the DC covariance matrix and assumed frequency-independent measurement error, which results in CDC = 2 DC , where σ DC is a standard deviation of DC data noise and δ ij is the Kronecker delta. The standard deviations (i.e., σ DC and σ GF ) are treated as hyper parameters and solved together with the model parameters (Bodin et al., 2012).

Prior Probabilities
We assume truncated uniform distributions for the prior probability of k, σ DC , and σ GF . We also assume the following limits: [k min , k max ) = [1, 51) for k, [ DCmin, DCmax] = [0.005, 0.090] for σ DC (unit in km/s), and [ GFmin, GFmax] = [0.02, 0.07] for σ GF . We tested several choices for these parameters to find that the resulting 7 of 15 velocity structures did not change significantly. We set the minimum limit of the layer depths to z min = 2.388 (water depth) and the maximum to z max = 15 km and use the Dirichlet partition prior with unit concentration parameters (Dosso et al., 2014). This setting corresponds to a non-informative prior: the prior probability remains constant no matter where the layer boundary is between z min and z max . We use the Gaussian distribution with a zero mean for the Vs anomalies. The Gaussian width (i.e., standard deviation σ δβ ) must reflect how reliable the reference model is. We set this parameter to 0.2 km/s. In summary, the joint prior can be expressed as follows: We confirmed that our inversion code implements the prior probability as intended by performing RJMCMC, forcing the likelihood to be zero ( Figure S12 in Supporting Information S1).

Probabilistic Sampling With Parallel Tempering
The RJMCMC algorithm aims to sample the posterior probability ( | ) through iteration. At each iteration, is proposed by either (a) adding a layer, (b) removing a layer, (c) moving a layer interface, (d) perturbing the S-wave velocity of a layer, or (e) perturbing the standard deviation of the data noise. One of the above-mentioned five procedures is randomly selected at each iteration to generate a new model. The proposed model is accepted at a probability α MHG , which is defined as the tempered Metropolis-Hastings-Green criterion (Green, 1995): is the probability that a transition from ( ) to is proposed, and |J| is the Jacobian compensating for a unit volume change in the model space. The exponent T(>1), which represents a temperature that loses the acceptance criterion, is a modification of the original Metropolis-Hastings-Green criterion. In the parallel tempering method (Geyer & Thompson, 1995;Sambridge, 2014), differently tempered Markov chains are run in parallel. At the end of each iteration, 10 pairs of chains are probabilistically allowed to swap the temperatures. Based on this swap, the random walk can undergo a long jump in the model space and efficiently converge to the global maximum.
The inversion involves 1,000,000 iterations, including the first 800,000 iterations of the burn-in period. In total, 100 Markov chains are run in parallel, 20 of which have a unit temperature and are used to evaluate posterior probabilities. We only save the models every 2,000 iterations to avoid artificial correlation between samples.

Results
The ensemble of model parameters sampled by the transdimensional inversion provides insights into the probable range of a 1D Vs structure beneath each station. Figure 5 shows the inversion results obtained at HDA06. The posterior marginal probability of Vs as a function of depth indicates a well-converged solution with a clearly defined peak at each depth. Other diagnostic information, such as the evolution of log-likelihood and acceptance ratio of proposals, is shown in Figure S13 and Table S2 in Supporting Information S1, respectively. According to the mode value at each 0.3 km depth (green line, Figure 5f), the velocity increases up to a depth of 4.8 km, with sharp, positive velocity contrasts at depths of 2.7 and 4.2 km. Although less clear, these discontinuities can be seen in the maximum a posteriori (MAP) estimate (purple line, Figure 5f). We conclude that these contrasts reflect different lithologies of sediments and refer to the layers as sedimentary units 1-3 (U1-3), from top to bottom.
Beneath this unit sequence, Vs abruptly drops to form a LVZ. The top of the LVZ is 0.1 km deeper than the depth at which the referenced Vp tomography model exhibits velocity inversion. Note that our prior Vs information already incorporates the velocity inversion that can be observed in the Vp model (black curve, Figure 5f). The AKUHARA ET AL.

10.1029/2022JB026298
8 of 15 inversion analysis requires the further reduction of Vs, suggesting a high Vp/Vs ratio in the LVZ: based on the assumption of a Vp of 3.4 km/s from the Vp tomography model, the Vp/Vs ratio corresponds to 2.8. However, this estimation likely overestimates Vp/Vs ratio. This is because the reference Vp tomography model has a coarser vertical resolution than Vs profiles obtained in this study, subject to smoothing constraints. Thus, we smoothed the Vs profile using a running window of 1.5 km depth to mimic the vertical resolution of seismic tomography ( Figure S14 in Supporting Information S1). The window length of 1.5 km was chosen by trial and error so that Vs profile exhibits a similar degree of smoothness to the reference Vp model. Even after this smoothing, the Vp/ Vs profile culminates at the LVZ with a maximum value of 2.5.
Inversion results from other stations show similar first-order features. Three layers (i.e., U1-3) are discernible immediately beneath the seafloor, and a LVZ can be detected beneath them, especially evident with mode estimations (Figure 6). An exception is HDA01 without a LVZ. This absence of LVZ beneath HAD01 could be artificial, considering that Vs profiles from the other stations consistently exhibit a LVZ. The overall consistency of the waveform to the other stations (Figure 4) suggests that a velocity model with a similar LVZ likely explains the HDA01 data well. Some undesirable features in the HDA01 data, perhaps originating from unmod eled components of actual velocity structure or high noise contribution, may distort the posterior probability estimation.

of 15
Since the LVZ is the center of interest, we exclude HDA01 results from the discussion for simplicity. Following the last paragraph, we calculate smoothed Vp/Vs profiles of each station. The resulting Vp/Vs profile shows a peak at the LVZ depth for most stations (Figure 7). The peak values from mode estimations are consistent among HDA02, 03, 04, 05, 06, 08, and 09, ranging from 2.5 to 2.7. Stations HDA07 and HDA10 show relatively lower Vp/Vs, 2.2 and 2.0, respectively, but the probability distribution of those stations has an elongated tail toward higher Vp/Vs. Thus, the Vp/Vs ratio of 2.5-2.7 may also be applicable to these two stations.
To quantify the depth of each lithological boundary, we searched for the depth of maximum velocity contrast within a given depth range. This search was performed for all 1D S-wave velocity structures sampled in the inversion. The aggregation of all results provides statistics for the lithological boundary depths. We set depth ranges for this search to 2.3-3.1 km for the boundary U1-U2, 3.1-5.5 km for U2-U3, 4.0-7.0 km for U3-LVZ, and 5.5-9.5 km for the bottom of the LVZ. The resulting median estimates are shown as background colors in Figures 6 and 7. In addition, 68% confidence intervals are shown in Figure 8b. Note that this error estimation tends to be biased toward magnifying uncertainties because the transdimensional inversion can produce ineffective (i.e., too thin) layers at random depths with a considerable velocity contrast. Hence, we chose to display the 68% confidence intervals in Figure 8b rather than the commonly used 95% intervals.
The above qualitative estimates of uncertainties confirmed the lateral variation in the depth of the top of the LVZ: the lithological boundary deepens on the southwestern side, whereas it becomes shallower on the northeastern side (Figure 8a). The depth offset is sharp: ∼1 km vertical offset within a distance of 0.5 km. The green vertical bars in Figure 8c show a 68% range of theoretical arrivals of the Ps converted phase from the top of the LVZ, The present study fixes a P-wave velocity structure at a single reference model and applies it to all stations, ignoring the presence of lateral heterogeneities and uncertainties in the reference model. However, such a fixed Vp could minorly bias Vs estimation because P-wave GFs (or receiver functions) have secondary sensitivity to Vp/ Vs ratios (e.g., Zhu & Kanamori, 2000). To quantify this effect, we solved Vp anomalies as well as Vs, where a Gaussian distribution with a standard deviation of 0.15 km/s is used as the prior probability for Vp anomaly. The results show that the main feature (i.e., the LVZ) does not change, irrespective of whether Vp is solved ( Figure  S15 in Supporting Information S1). The posterior probability of Vp remains nearly identical to the prior probability below a depth of 4 km, suggesting that the data set is only sensitive to the shallow part of the Vp structure. The longer time window of GFs could constrain Vp/Vs ratios of the LVZ, but unfortunately, GFs do not show good consistency for phases arriving later than 4 s (see Figure 4a).
Another concern is overfitting. The transdimensional inversion could unnecessarily add many thin layers to cause overinterpretation of input data. In theory, this issue can be avoided by the adopted transdimensional inversion scheme but could occur with an inappropriate parameterization made for the likelihood, for example. To see whether the obtained LVZ is robust, we enforced a smaller number of layers by setting k max to 21. Still, we observe an evident LVZ ( Figure S16 in Supporting Information S1). As another test case, we conducted a fixed-dimensional inversion by fixing k at 20. The other parameters, including layer depths, are allowed to vary freely. We found that this fixed-dimensional setting fails to reach a well-converged solution, highlighting the efficient model search by the transdimensional algorithm ( Figure S17 in Supporting Information S1). This kind of advantage in the transdimensional scheme has not been discussed elsewhere, to the best knowledge, but should be investigated more in the future.

Geological Interpretation
The inversion results present a remarkable low-velocity, high Vp/Vs feature with a velocity inversion. Typically, marine sediments undergo a monotonic increase in Vs with increasing depth because of compaction 12 of 15 (Hamilton, 1979). The velocity inversion observed in this study is unexpected. A plausible cause for the observed velocity inversion is high pore fluid pressure. Based on theory and experiments, it is known that high pore fluid pressure increases the Vp/Vs ratio of marine sediments (Dvorkin et al., 1999;Prasad, 2002), which agrees with our results. Therefore, we interpret that the LVZ represents a fluid reservoir. Aligned cracks could also explain the high Vp/Vs ratio even in the absence of fluid through anisotropic effects (X. Q. Wang et al., 2012). However, we find that numerical modeling based on a scattering theory with penny-shaped parallel cracks (Hudson, 1981) fails to explain such high Vp/Vs ratios (2.5-2.6), at least within the reasonable range of crack density (<0.1; Crampin & Leary, 1993). It only predicts a Vp/ Vs ratio up to 2.0. For this modeling, we assume an isotropic host rock with a Vp of 3.6 km/s, Vs of 2.0 km/s, and density of 2.3 g/cm 3 . Those values are extracted from Unit 3.
Sustaining the overpressure condition within the fluid reservoir will require a relatively impermeable structure above. Laboratory measurements on terrigenous sediments from deep-sea drilling have shown that the porosities gradually decrease with depth, from ∼70% at the sea bottom to ∼20% at a burial depth of 1.5 km (Kominz et al., 2011). It has also been reported that porosity changes from 70% to 20% for mudrocks correspond to 3-4 orders of magnitude decrease in permeability (Neuzil, 1994). Thus, we speculate that the bottom of Unit 3, with a burial depth of ∼2.6-3.9 km, undergoes more severe porosity loss and can impede fluid to permeate shallower layers. This permeability barrier could trap abundant fluid below, leading to the formation of the fluid reservoir.
Considering the shallow subduction depth (∼10 km), subducted sediment along with the Phillippine Sea plate is likely a fluid source, which can release fluid via mechanical compaction or dehydration (Saffer & Tobin, 2011). The occurrence of slow earthquakes may reflect fluid-rich conditions near the subducting plate interface: lines of evidence require high pore fluid pressure for the genesis of slow earthquakes (Behr & Bürgmann, 2021 and references therein). Since this possible fluid source is vertically separated from the LVZ, permeable structures such as faults or fractures will be required to effectively convey fluids from the subducted sediments to the LVZ, as discussed in the next paragraph. Such permeable structures may not penetrate Unit 3. After reaching the bottom of Unit 3, fluids might diffuse laterally in accordance with permeability anisotropy due to sediment stratification.
The presence of faults in the overriding prism seems natural for this region with the subducted KPR. Analog and numerical experiments have demonstrated that many back-thrusts occur on the leading flank of the seamount (Dominguez et al., 1998;Sun et al., 2020). A recent compilation of seismic reflection surveys in the Hyuganada has identified several NNW-SSE trending thrust faults northeast to the array (Figure 9; Headquarters for Earthquake Research Promotion, 2020). At ∼2 Ma, before the last change in the convergence direction, the KPR was located east of the present position (Mahony et al., 2011). The subsequent oblique subduction involves right-lateral motion along the trench, potentially inducing the northeast-dipping back-thrust near the array (Figure 10a). Notably, this fault trend is roughly parallel to the sharp offset in the LVZ depth we observed. The sharp offset could imply the presence of a blind back-thrust fault beneath it. Cumulative deformation along the fault might be responsible for the sharp offset. If existing, such a fault will act as a fluid conduit (Figure 10a).
We acknowledge that our data set poses only weak constraints on the geological process behind the LVZ and thus does not exclude other possibilities. For alternative hypothesis, the LVZ could represent ascending, overpressured material, such as a mud diapir (e.g., Brown, 1990), about to pierce into Unit 3 (Figure 10b). The head of ascending body would selectively intrude into Unit 3 along a mechanically weakened fabric parallel to the KPR, which leads to the NNW-SSE trending depth offset. Similarly oriented faults nearby the array (Figure 9) support the presence of such a weak fabric.
Another alternative is intense fracturing (Figure 10c). While our conservative modeling using parallel, isolated cracks excludes the possibility of the absence of fluid, a more severely damaged medium by fractures and consequent seismic anisotropy (e.g., Schoenberg & Douma, 1988) could explain the high Vp/Vs ratio without fluid. One advantage of this scenario is that there is no requirement for impermeable layering at the base of Unit 3. For the hypothesis by a fluid reservoir, it is still unclear whether such an impermeable layer can remain intact under severe tectonic deformation caused by ridge subduction. A drawback, on the other hand, is a lack of reasonable explanation why the LVZ develops almost horizontally over a wide area. In order to realize the high Vp/Vs ratio, fractures are likely to be aligned with their plane extending vertically (Ding et al., 2019), which seems somewhat unsuited to forming the horizontally flat LVZ. Further investigation in combination with active-source seismic surveys can illuminate the cause of the LVZ, but this is left for our future work.
This study has identified that the LVZ extends laterally, at least to the array size (∼4 km). The Vp gradient profile shown in Figure 1d suggests that the LVZ extends ∼60 km laterally beyond the aperture of the OBS array. Moreover, an independent seismic refraction profile in Hyuga-nada has obtained a comparable low-velocity feature within the accretionary prism, ∼50-100 km south of the array (Nishizawa et al., 2009), possibly suggesting that similar fluid reservoirs are widely distributed in this region. Pursuing its spatial extent will be important for better understanding the cause of the LVZ and hydrological processes of Hyuga-nada in association with the KPR.

Conclusions
In this study, the Vs structure in the Hyuga-nada accretionary prism was constrained using a passive seismic array. The Vs structure exhibits a LVZ beneath stratified sedimentary units (U1-3). Based on the reduced Vs and high Vp/Vs ratio, we conclude that the LVZ reflects a fluid reservoir with high pore fluid pressure sustained by the impermeable layering above. The significant depth offset of the top of the LVZ, extending over ∼4 km of the array aperture, possibly suggests the presence of a blind thrust fault. Such faults generated by the subduction of the KPR may act as fluid pathways and contribute to the reservoir. However, we do not exclude other possibilities: the LVZ may reflect a mud diapir or intense fracturing, for example. 14 of 15 The results of this study demonstrate the potential of passive seismic source analyses to acquire a high-resolution structure of Vs, leading to gaining new constraints on fluid processes in the accretionary system. A limitation is its narrow resolvable range laterally, which may hamper interpreting resultant Vs structures conclusively. Joint interpretation with active seismic source surveys will remedy this drawback. Nowadays, a number of seismic survey data have been obtained in subduction zones worldwide. Additional passive seismic experiments like this study will help understand physical properties and hydrological features in the accretionary prism.