Influence of a Subducted Oceanic Ridge on the Distribution of Shallow VLFEs in the Nankai Trough as Revealed by Moment Tensor Inversion and Cluster Analysis

The eastern Nankai Trough is a unique site where many shallow, very low frequency earthquakes (sVLFEs) are recorded by nearby broadband ocean bottom seismometers. Here, we estimated the locations and seismic moment tensors (MTs) of sVLFEs based on the low‐frequency (<0.07 Hz) components of the records. Although some sVLFEs exhibited long‐duration signals (>100 s), indicating a degree of source complexity, the MT inversions were limited to events with impulsive and short durations (20–30 s). Nevertheless, cluster analysis confirmed that the MTs of impulsive sVLFEs reasonably represented the overall deformation for an event group including many longer‐duration events. The distribution of MTs indicates that deformations associated with sVLFEs are influenced by a subducted oceanic ridge in this region, which produces an along‐the‐trough variation in a strain accumulation/release pattern and probably controls the spatial patterns of tsunami generation.


Introduction
Shallow, very low frequency earthquakes (sVLFEs) are anomalous earthquakes of a long duration and signals depleted at high frequencies relative to most regular earthquakes of the same magnitude. They have been detected in some shallow (depth <~15 km) subduction zones (reviewed in Obara and Kato, 2016). The eastern Nankai Trough, where the Philippine Sea plate subducts beneath the Eurasian plate, is a unique site where sVLFEs have been recorded just above their sources by broadband ocean bottom seismometers (BBOBSs). The first such successful close-in observation was made by Sugioka et al. (2012) using an array of three temporarily deployed BBOBSs (Figure 1a). Later on, more observations have been made over a broader region by a cabled ocean floor network DONET1 (Kaneda et al., 2015;Figure 1a). Based on the results from full-waveform moment tensor (MT) inversions of low-frequency components (<0.125 Hz), sVLFEs are considered to be (i) generated by low-angle thrust faulting in the décollement zone along the plate boundary and (ii) laterally spread from beneath the outer ridge toward the deformation front (Nakano et al., 2016(Nakano et al., , 2018Sugioka et al., 2012). Alternatively, based on the analysis of high-frequency components (>2 Hz) of sVLFEs recorded by DONET1, Toh et al. (2018) suggested a possible along-the-strike variation in their distribution related to the subduction of the Paleo-Zenisu ridge (Park et al., 2003;Figure 1a), which is located on the northeastern side of DONET1. However, the high-frequency components of sVLFEs are strongly scattered by heterogeneous structures; thus, estimated source locations have large uncertainties. Consequently, the key observations were limited to locations estimated for only several sVLFEs, whose onset of direct P arrivals was recognizable, out of more than a hundred events. Nevertheless, revealing the sVLFE distribution with respect to the subducted ridge, based on more substantial amount of data, should help understand its influence on the stress field and its strain accumulation/release pattern, which can also link to spatial tsunami generation patterns. Furthermore, it can tell us about the physical conditions that promote the sVLFEs.
Here, we estimate the locations and focal mechanisms of sVLFEs using the low-frequency components (<0.0625 Hz) of DONET1 BBOBS records. An overall distribution was estimated by combining results from a full-waveform MT inversion applied to a limited number of sVLFEs and a cluster analysis based on amplitude patterns measured for all the detected events. Finally, we discuss the influence of the subducted ridge to the sVLFE distribution.

Data
We used three-component ground velocity records from 18 BBOBSs (Guralp CMG-3T) that were equipped in DONET1. Out of the 20 stations, Stations C10 and C11 were not used due to their high noise levels. We used the records for sVLFE activity that continued for 1-3 weeks each year immediately beneath the network during 2015, 2016, and 2017, while the previous study was limited to the data from 2015 (Toh et al., 2018). Using the classification of the local seismicity beneath DONET1 (supporting information Table S1; To et al., 2015), we visually picked 335 VLFEs whose signals were observed by at least six stations on the vertical component filtered between 16 and 32 s.

Event Selection
While some sVLFEs showed somewhat impulsive signals in the low-frequency range (e.g., Figure 2a), others exhibited complex waveforms with long durations (e.g., Figure 2b). These complexities were likely produced by their long source time functions (Sugioka et al., 2012) and the spatial complexities of the source, both of which require a spatiotemporal finite-source method to resolve. In this study, we focused on the impulsive sVLFE signals and applied a simpler spatiotemporal point source model. Also, because sVLFEs tend to occur in spatial clusters (e.g., Toh et al., 2018), the MT inversion was applied only to those apparently spatially clustered events and not applied to sparsely located ones. Merits of imposing this limitation in the analysis will be explained in section 6.2. Conversely, an amplitude pattern was measured for all the detected events. The pattern can be measured for any event yet reflects the location and focal mechanism of the event and can be used as their proxy. Then, a cluster analysis was performed to systematically evaluate the amplitude patterns and identify groups of closely located events. The amplitude pattern was obtained by measuring the maximum amplitude from the envelope of the vertical velocity seismograms, which were obtained after removing the instrument response and bandpass filtering at 16-32 s. This was measured at a time up to 65 s after the earliest signal peak in the entire network. The maximum amplitude was measured when the signal was visually observed or when it was below a certain background noise level (<5.0 × 10 −8 m/s; i.e., seeing no signal is also an important information) and not when a large noise obscured the signal. Figure 1b presents the amplitude patterns for all 335 events. Similar measurements for teleseismic events, listed in Table S2, were nearly flat ( Figure 1c). Hence, the significant variations in the observed patterns for the sVLFEs were not due to site responses. Many amplitude patterns displayed a sharp peak at a single station, which indicates that the station was located well within the near field of the sVLFE source.
Clustering was conducted using HDBSCAN (Campello et al., 2013;McInnes et al., 2017), which is a density-based hierarchical clustering method that allows for the extraction of clusters of different density thresholds. The pairwise distance (d AB ) (or dissimilarity measure) between any two events (A and B) is defined by Equations 1-5: where a A and a B are the amplitude pattern vectors measured for events A and B, respectively, and a min is the minimum amplitude detected among all the events. Therefore, x A and x B consist only of nonnegative elements, and small amplitude measurements (i.e., noise) are converted to values close to 0. Both c A and c B are cosine similarities between x A and x B . The number of elements (i.e., stations) used to evaluate the cosine similarities was set to 6. The i A and i B are sets of six stations in descending order of amplitude observed for events A and B, respectively, among the stations commonly available between events A and B. After cluster analysis, 197 of the 335 sVLFEs were assigned to 26 groups (G00 to G25) that consisted of 2 to 35 events. Figure 3 plots the amplitude patterns for each group.
The application of the MT inversion was limited to the 197 grouped events. Furthermore, it was limited to 134 events with relatively impulsive waveforms judged by visual inspection. Our strategy is to take advantage of the nature of sVLFEs, whereby events of simple (e.g., Figure 2a) and complex (e.g., Figure 2b) waveforms were often mixed within the same spatial cluster. For example, the two events in Figure 2 were assigned to G11 and showed a similar amplitude pattern (Figure 3). If at least one event per group is resolved by the MT inversion, we will obtain the overall distribution of the spatially clustered sVLFEs. . Instrument responses were removed from the observations. Both were filtered at 16-32 s. Amplitudes were normalized by the maximum observed amplitude per station. The traces of malfunctioning instruments were not plotted. The event location was estimated to be right below Station D13 (yellow beachball in Figure 4a). Horizontal components were rotated to the tangential and radial components based on the estimated location. (b) Observed waveforms for an sVLFE on 7 April 2016. The horizontal components were rotated by assuming the same location as in (a).

Geophysical Research Letters
where the weighting factor, w i , is set proportional to the root square of the distance between the source and the station. We weighted the vertical components five times more than the horizontal components because they are known to be less noisy (Suetsugu & Shiobara, 2014), yet their signal amplitudes are much smaller than those of the horizontal components (e.g., Figure 2). Both the observed and synthetic waveforms were bandpass filtered either between 16 and 32 s or between 20 and 50 s, depending on the signal duration and signal-to-noise ratio per event. The location and time were grid searched at intervals of 0.01°i n the horizontal direction, 1 km in depth between 5 and 11 km below sea level (bsl), and 0.5 s in time. For each grid, the MT was estimated using the code TDMT_INV (Dreger, 2003) with some modifications made to the weighting scheme, as previously described. The code conducts a full-waveform linearized inversion in the time domain. Green's functions were calculated using the software developed by Herrmann (2013), which is based on the frequency wavenumber integration. We used a 1-D P wave velocity (Vp) model (shown in Figure S5a of Text S1) that is constructed by horizontally averaging a 2-D model obtained by an active source survey (Nakanishi et al., 2008). The S wave velocity and density profiles were constructed based on empirical relationships to Vp (Brocher, 2005). The variations in the station elevation (1.9-3.8 km bsl) were considered by preparing Green's functions for three different station elevations (i.e., water layer of 2.00, 2.40, and 3.65 km bsl) and by using the one closest to the actual station elevation.

Results: MT Inversion
We obtained 65 sVLFEs (Figure 4a) that satisfied the following two criteria: (i) a VR of >65% and (ii) the estimated centroid source time was less than 6 s before the earliest peak arrival observed in the entire network (explained in Figure S1b). The colors of the beachballs in Figure 4 represent the sVLFEs groups. The criteria were chosen so that events for the same group were closely located. By lowering the criteria, outliers began to appear in some of the groups ( Figure S1). The standard errors of source locations estimated by the bootstrap method for three arbitrarily selected sVLFEs were 2 and 1 km in the horizontal and vertical directions, respectively ( Figure S2). The obtained moment magnitudes ranged from 2.4 to 3.9 ( Figure S3). Almost all the VLFEs could be reasonably explained by the low-angle thrust faulting, as previously reported. Some comparisons of the synthetic and observed waveforms are provided in Figure 2a and in Figure S8 of Text S2.

Geophysical Research Letters
TOH ET AL. Figures 4b and 4c plot the sVLFEs along the cross sections for the northeastern and southwestern areas. We note that events were plotted at a depth 1.5 km shallower (indicated as "adjusted depth") than originally obtained, such that the depth distribution appears consistent with the results of Sugioka et al. (2012). A possible cause of this apparent systematic shift in the estimated depth is discussed in Text S1. Figure 4d shows the total number of events per group along with beachballs that represent the number of resolved events per group. No events were resolved for G03, G19, and G22. The events in G03 were unresolved because they occurred close in time with other events. Nevertheless, the sharp peak at Station B06 in the amplitude patterns (Figure 3) indicates that they occurred close to this station.

Geophysical Research Letters
Events in G19 and G22 likely occurred outside of DONET1 because their amplitude patterns exhibited no clear peak. At least one event was resolved for each of the other groups. Thus, the distribution of the clustered events is represented by Figures 4a-4c. 6. Discussion and Conclusions

sVLFE Less Active Region on the Northeastern Side
Figures 4a-4c suggest that the sVLFEs occurred at different distances from the deformation front between the northeastern and southwestern sides. Hereafter, the region with no large sVLFE on the northeastern side of DONET1 is referred to as the sVLFE, or seismically, less active region (SLA region; gray ellipse in Figure 4a). The SLA region could also be an artifact caused by the heterogeneous station distribution. Nevertheless, we verified that it is the actual feature of the sVLFE distribution as follows.
First, we demonstrated that the observed and synthetic amplitude patterns were reasonably comparable. Figure 5b shows the synthetic amplitude patterns for the 65 resolved events, as obtained by constructing synthetic waveforms based on the MT inversion results. Figure 5a presents the observed amplitude patterns of the corresponding events. Comparing Figure 5a with 5b reveals that the observations were reproduced roughly well by the synthetics.
Second, we showed that the lack of sVLFEs in the SLA region was not due to the elimination of such events during data processing. To demonstrate this, synthetic events with low-angle thrust faulting mechanism were placed in the SLA region (Figure 5e). The predicted amplitude patterns for the synthetic events are plotted in Figure 5d. Meanwhile, observed amplitude patterns for all the detected events whose peak amplitude was observed at Station B06, including those for which the MT inversion was never applied, are plotted in Figure 5c. The amplitude at B06 is anomalously large compared to other stations in the observations (>5 times the average; Figure 5c), whereas it is not that large in the synthetics (<4 times the average; Figure 5d). There is one exception in the observations, whereby the peak at B06 was only twice the average; however, we checked and confirmed that the signal of this event was contaminated by surface waves of a teleseimic event.
The synthetic test showed that the synthetic events located in the SLA region produced amplitude patterns that were never observed and, thus, suggests that no events were detected in the region.
Finally, we showed that the absence of sVLFEs in the SLA region was not due to the detection ability of DONET1. Figure 5f plots the unnormalized amplitude patterns for the synthetic events along with the background noise. Synthetic events were assigned an Mw of 3.0, which is comparable to detected events ( Figure S3). The amplitude patterns for the background noise were first measured for the 335 time windows that start at 200 s before the earliest arrivals of the detected sVLFEs and last 65 s. Out of these 335 amplitude patterns, 153 patterns of relatively quiet time with a maximum amplitude of <1.5 × 10 −8 m/s at all the available stations are plotted in Figure 5e. Figure 5e suggests that the sVLFEs with Mw ≥ 3 should be detected by observing their signals at more than six stations, which was our event detection criterion. Thus, the synthetic test suggests no sVLFEs with similar or larger magnitudes occurred there. Comparisons between the noise and synthetic signals are further described in Text S2, and comparisons with other studies are provided in Text S3.

Why Perform Cluster Analysis?
Compared to regular earthquakes, less is known yet about the source process of sVLFEs. Here, we selected events with simple waveform and evaluated the modeling by VR, which does not distinguish between (1) noise in the data and (2) possible flaws in the theory, such as the point source or structural assumptions. Thus, groups of events that interest us, but whose source process deviates from the assumptions, may be eliminated during the data processing. The cluster analysis based on amplitude patterns was introduced to identify such event groups and to ensure that the obtained results represent the overall distribution. Fortunately, the groups eliminated in this study were only those likely located far from the network (i.e., G19 and G22). Conversely, some interesting observations were elucidated by the analysis. For example, an event in G23 shows a focal mechanism, unlike the others and the cause of which is unknown yet (Figure 4d). Also, Figures 5a and 5b suggest that amplitude patterns for some groups with their peak and a second peak at D14 and D15 (e.g., G17 and G18) are systemically different between observations and synthetics. Since these are events with relatively impulsive signals, the discrepancy is likely caused by factors 10.1029/2020GL087244

Geophysical Research Letters
other than the source duration, such as the finite fault size (Ito & Obara, 2006b) or structures. Future work should, therefore, invoke more complex source models to resolve the sourcing process for a larger number of sVLFEs.

Subducted Oceanic Ridge Influence on sVLFEs Distribution
The sVLFEs are considered to occur at fault zones that have transitional frictional properties between stick slip and steady aseismic slip. Further, high pore fluid pressure can be one of the causes for such transitional state to arise (Saffer & Wallace, 2015;Schwartz & Rokosky, 2007). Here, we discuss that the subducted ridge likely produces preferential conditions for sVLFEs occurrence on its downdip side and thus reasonably explains the obtained distribution.
On the southwestern side of DONET1, sVLFEs have been known to distribute from below the outer ridge toward the deformation front (Nakano et al., 2018;Sugioka et al., 2012), and our results are consistent (Figures 4a and 4c). The décollement zone along the plate boundary is considered as an updip pathway for fluid generated in deeper parts of the subduction zone (Brown et al., 1994), and it has been suggested that the sVLFE locations coincide with an area of high pore fluid pressure along the décollement zone (Kitajima & Saffer, 2012). Also, the sVLFEs were interpreted to indicate possible tsunamigenic slip along the plate boundary, extending toward near the deformation front (Sugioka et al., 2012).
Our results suggest that the story is different on the northeastern side of DONET1 due to the ridge subduction, as also described by Toh et al. (2018). The sVLFEs are concentrated downdip of the peak of the subducted ridge (Figure 4b). We speculate that the pore pressure should be high in this area, because the fluid from the downdip side would stay there, rather than to pass over the ridge of low permiability to follow the décollement zone trenchward. Also, larger shear stress would be exerted ahead (i.e., downdip) of the topographic high rather than behind it (Dominguez et al., 1998). The tsunamigenic potential may be high along imbricate thrust faults, located~25 km landward of the deformation front, that extend upward ( Figure S4b, blue line) from the downdip flank of the subducted ridge where the strain accumulation/release appears the most active. Finally, our observations agree with the conceptual model that subducting seamounts would develop a fracture network in the upper plate to host small earthquakes and aseismic creep, rather than to nucleate large earthquakes (Wang & Bilek, 2011). Such a fracture network could be the cause of the transitional frictional properties of the fault.

Data Availability Statement
DONET1 waveform data are available via the website of the National Research Institute for Earth Science and Disaster Resilience, Japan (http://www.hinet.bosai.go.jp/). Data Set S1 is available online (at https:// zenodo.org/record/3928970#.Xv64IfLgoUB).