Locating tectonic tremors with uncertainty estimates: Time- and amplitude-difference optimization, wave propagation-based quality control, and Bayesian inversion

Summary The accurate location of tectonic tremors helps improve understanding of their underlying physical processes. However, current location methods often do not statistically evaluate uncertainties to a satisfactory degree and do not account for potential biases due to subsurface structures not included in the model. To address these issues, we propose a novel three-step process for locating tectonic tremors. First, the measured time- and amplitude differences between station pairs are optimized to obtain station-specific relative time and amplitude measurements with uncertainty estimates. Second, the time– and amplitude–distance relationships in the optimized data are used to roughly estimate the propagation speed (i.e., shear wave velocity) and attenuation strength. Linear regression is applied to each event, and the resulting velocity and attenuation strength are used for quality control. Finally, the tremor location problem is formulated within a Bayesian framework where the model parameters include the source locations, local site delay/amplification factors, shear wave velocity, and attenuation strength. The Markov chain Monte Carlo algorithm is used to sample the posterior probability and is augmented by a parallel tempering scheme for an efficient global search. We tested the proposed method on ocean-bottom data indicating an intense episode of tectonic tremors in Kumano-nada within the Nankai Trough subduction zone. The results show that typical location errors (1 𝜎𝜎 ) are ~1–2 km horizontally and <5 km vertically. A series of experiments with different inversion settings reveals that adopting amplitude data and site correction factors help reduce random error and systematic bias, respectively. Probabilistic sampling allows us to spatially map the probability of a tremor occurring at a given location. The probability map is used to identify lineaments of tremor sources


Introduction
Tectonic tremors, considered as a swarm of low-frequency earthquakes, constitute a broad spectrum of slow earthquakes together with very low-frequency earthquakes and slow-slip events. They were first discovered in southwestern Japan (Obara, 2002) and have since been identified at subducting plate interfaces worldwide (Araki et al., 2017;Brown et al., 2005;Nishikawa et al., 2019;Payero et al., 2008;Plata-Martinez et al., 2021;Rogers, 2003;Todd et al., 2018;Yamashita et al., 2015). Slow earthquakes, including tectonic tremors, release seismic energy over a long time considering their magnitudes, which indicates that they may be governed by different physical processes than regular earthquakes (Ide et al., 2007). Owing to their proximity to the rupture areas of megathrust earthquakes, slow earthquakes have drawn significant attention for their potential to deepen our understanding of future devastating earthquakes (Obara & Kato, 2016).
The accurate location of tectonic tremors is vital to understanding the slip behavior of plate interfaces. The spatiotemporal evolution of tectonic tremors has several unique but ubiquitous characteristics. First, tremors occur episodically, with their epicenters migrating parallel to the subduction margin, which indicates the simultaneous occurrence of slow-slip events. Second, tremors occasionally back-propagate against their main front at distinctly high speeds, known as rapid tremor reversal (e.g., Houston et al., 2011). Third, streaks of tremors in the dip direction of the subducting plate have been observed (e.g., Ghosh et al., 2010). These spatiotemporal patterns of tremors can constrain the frictional properties of the plate interface (Rubin, 2011), underlying physical processes (Cruz-Atienza et al., 2018), and structural factors that cause tremors (Ide, 2010).
The signals of tectonic tremors emerge without a clear phase onset, which makes locating them using the same methods as for regular earthquakes impractical. A common approach is the envelope correlation method (e.g., Mizuno & Ide, 2019;Obara, 2002), which cross-correlates enveloped seismograms between pairs of stations and assumes that the resulting time lag represents a difference in S-wave travel time. Optimization methods can then be applied to determine the source locations that best explain the measured arrival time differences. Another approach is to use the amplitude (e.g., Husker et al., 2012;Ogiso & Tamaribuchi, 2022), although such techniques are more widely used for locating volcanic tremors rather than tectonic tremors. Because seismic waves lose energy This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International during propagation, the spatial pattern of amplitudes can provide clues about source locations. However, this approach requires knowledge of attenuation structures and local site amplification, which typically necessitates additional analysis. Some studies have used a joint approach that combines both time-and amplitude-based methods, where the different datasets are often weighted subjectively (Maeda & Obara, 2009).
Despite the importance of investigating the source locations of tectonic tremors, many studies have not formally estimated the uncertainties associated with these locations, which raises the risk of misinterpreting results. Accurately estimating the uncertainties of tremor locations requires considering the statistics of the input measurements in the data domain (i.e., time and amplitude domains) and then converting them into the spatial domain by forward calculation. Uncertainties in the structure model used for the forward calculation must also be considered to prevent systematic biases. Such uncertainties in structures would be severe for offshore studies targeting shallow tectonic tremors (e.g., Yamashita et al., 2015). Typically, the seafloor is covered with unconsolidated sediments.
Such near-surface structures amplify the amplitude and delay the arrival of seismic waves, and the degree of this effect varies according to the geographic location.
To address the above issues, we propose a three-step method for locating tectonic tremors and estimating their uncertainty, which we applied to real tremor data obtained at Kumano-nada in the Nankai Trough subduction zone as a demonstration.

Data
We collected data from a seismic network at Kumano-nada in the Nankai Trough subduction zone, where the Philippine Sea plate subducts beneath the fore-arc margin. As shown in Fig. 1, the network comprises16 permanent cabled stations from the Dense Ocean Network for Earthquake and Tsunamis (DONET) (Kaneda et al., 2015;Kawaguchi et al., 2015) and 15  Intense episodes of slow earthquakes, including tectonic tremors and very lowfrequency earthquakes, repeatedly occur in this region at intervals of ~5 years (e.g., Takemura et al., 2022). The latest and most intensive episode began on December 6, 2020, and persisted for approximately 2 months (Ogiso & Tamaribuchi, 2022) within the observation period of the OBSs. We collected data from a 20-day period in the middle of this episode from January 1 to January 21, 2021, during which tremors are known to have occurred near the center of the seismic network (Ogiso & Tamaribuchi, 2022).
Note that our study focuses on accurately locating tectonic tremors rather than attempting to locate them exhaustively. This is why we only used 20 days of data, rather than the ~2 months of data covering the entire tremor episode, to reduce computation time.
For the same reason, we only briefly describe the tremor detection process below.
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International To detect the tectonic tremors, we preprocessed continuous seismic waveform data as follows. First, 300-s time windows were successively extracted from the continuous data with 50% overlap. The extracted time series were de-trended, tapered, bandpassfiltered from 1-10 Hz, and converted to envelopes via the Hilbert transform. We then smoothed the resulting envelopes with a 6-s triangular filter and merged the two horizontal components by using the root sum squared method. We did not use the vertical component because shear waves dominate the seismic records of tectonic tremors.
We detected tremors by inter-station cross-correlation. We used the 300-s envelopes to calculate cross-correlation functions between all possible station pairs every 150 s. For each station pair, all obtained correlation values with lag times from -150 s to 150 s were saved over the whole 20-day period to generate a histogram specific to the station pair.
We defined tremor detection as when a correlation value exceeded the 98th percentile of the histogram for at least 300 station pairs. Of the 11,520 time windows, 4324 time windows met these criteria.

Method
Our proposed method has three steps.
Step 1 is to optimize measurements from station pairs such as the arrival time difference and logarithmic amplitude ratio, which outputs the relative arrival time and logarithmic amplitude ratio at each station along with their respective uncertainties. These uncertainties can be incorporated in the final inversion stage to acquire uncertainties in the spatial domain.
Step 2 is to extract the firstorder features of wave propagation from the optimized station-specific data: the propagation speed and attenuation strength. These features are then used as quality control factors to retain good-quality data.
Step 3 is to invert the station-specific data and their uncertainties for hypocenters by using the Markov chain Monte Carlo (MCMC) algorithm in a Bayesian framework. To address biases from unknown structures, we jointly solve multiple hypocenters and include structural parameters and the associated correction factors in the model parameters.

Step1: Optimization of arrival time and amplitude differences
The unclear phase onset makes direct measurements of the arrival times of tectonic tremors a challenge. A widely used alternative approach is to use cross-correlation to measure the arrival time difference between station pairs (e.g., Obara, 2002): where ( ) is an envelope waveform recorded at the th station and Δ is the arrival This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International time difference between the th and th stations. This approach only works when the two waveforms are sufficiently similar. If the waveforms differ (e.g., due to different propagation paths), the measured arrival time difference can deviate from the true value.
In addition, a high level of noise can easily pose artificial peaks in the cross-correlation functions. Once the arrival time difference is obtained, the amplitude ratio between the two envelopes is defined as follows: This definition corresponds to the maximum likelihood estimation (MLE) of the amplitude ratio between two similar waveforms (Appendix A). The numerator has already been calculated to find the maximum of the cross-correlation function in Equation (1), so it does not require additional computation. Other definitions than Equation (2) may be used for the amplitude ratio, such as the squared sum (Maeda & Obara, 2009) or median value (Li et al., 2022). The obtained amplitude ratios are converted to amplitude differences by taking the logarithm so that they can be treated mathematically in the same manner as the arrival time differences: The above process yields ( − 1)/2 pairs of measurements, where is the number of stations. Individual pair measurements are dependent on other pairs (i.e., Δ ∼ Δ − Δ ). In other words, the ( − 1)/2 measurements inherently include redundancy. We may optimize this redundancy by solving a linear system (VanDecar and Crosson, 1990): where 1 ⋯ denote the relative arrival time at each station. A regularization condition is added to the bottom row that imposes a zero-sum requirement on the relative arrival times. This system has the following analytical solution (VanDecar & Crosson, 1990): This optimization reduces the redundant measurements of ( − 1)/2 station pairs to individual station-specific measurements. The original redundancy reflects This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International the measurement errors. The corresponding standard deviation can be calculated as follows (VanDecar & Crosson, 1990): The later Bayesian inversion in Step 3 takes the above station-specific relative arrival times ( ) as input data, and their standard deviation ( ) are used for calculating the likelihood.
Similar equations hold for logarithmic amplitudes: where is a relative logarithmic amplitude at the th station, and is the corresponding standard deviation. Such poor-quality data, even if present at only a few stations, can distort the optimized solution significantly because the optimized solutions given by the arithmetic mean (i.e., Equations (5) and (7)) are not robust against outliers. This sensitivity to poor-quality data requires an automatic and objective process to reject ill-optimized results.
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International

Step2: Wave propagation-based quality control
The optimization in Step 1 is useful for capturing seismic wave propagation intuitively. In cases where the optimization is successful and not affected by outliers, the relative arrival times and amplitudes exhibit a concentrated pattern when viewed on a map where the center approximates the epicenter, as shown in Fig. 3. We can use this pattern to obtain time-distance and amplitude-distance relationships, which in turn can be used to roughly quantify the propagation speed (i.e., S-wave velocity ) or attenuation strength (i.e., quality factor ), respectively.
For a uniform velocity structure throughout the medium, the arrival time is proportional to the propagation distance : This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International Thus, can be estimated from the slope of the time-distance plot (Fig. 3c).  (6) and (8).
(c, f) Distance plots of the optimized station-specific measurements. The error bar denotes the standard deviation. The blue dashed line represents a regression line.
The amplitude of a body wave at a propagation distance is described as where 0 is the source amplitude, is the representative frequency, and is the quality factor. Taking the logarithm of Equation (10) After correcting the geometrical spreading (i.e., adding ln to Equation (12)), the logarithmic amplitude becomes proportional to the distance. Therefore, we can determine the attenuation strength from the slope of the amplitude-distance plot (Fig. 3f).
We propose using the estimated and (or if and are fixed) as quality control factors to select good-quality events. The estimated and are representative of a broad region where source-receiver paths pass through. Because This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International tectonic tremors always occur in a narrow depth range, all ray paths most likely propagate through similar depths. Considering that subsurface properties vary less laterally than vertically, the estimated and values from different events should fall into a narrow and physically reasonable range. Hence, events with outlier and values can be attributed to ill-optimized datasets or events far isolated from target tremors, such as teleseismic events.
In practice, the propagation distance is not known before the hypocenter is determined. Instead, we may assume that the source is located beneath the station with the maximum relative amplitude. In this study, we assume that the focal depth is 7 km below sea level considering the specific tectonics of the study area. The resulting distance plots from the assumed source location are then linearly regressed by the least squares method. The resultant and values from different events are shown in Fig. 4. Based on the scatter plot, we selected events with =2.0-4.0 km/s and =0.015-0.030 as acceptable. These ranges are comparable to those previously estimated for the study area (Akuhara et al., 2020;Yabe et al., 2021), and they correspond to of 130-520 if a dominant frequency of 5 Hz is assumed. Under this criteria, 580 of the 4324 events were retained.

Step 3: Bayesian inversion
In Step 3, we adopt a Bayesian interface to invert the relative arrival times and logarithmic amplitudes jointly for the hypocenters ( , , ; = 1, ⋯ , ) , delay factor for each station ( ; = 1, ⋯ , ), amplification factor for each station ( ; = 1, ⋯ ), S-wave velocity ( ), and quality factor ( ). Here, and represent the numbers of stations and events, respectively. The delay and amplification factors are used to account for the local effects caused by seafloor sediment beneath the stations. We assumed uniform structures for the S-wave velocity and attenuation for simplicity. These model parameters are denoted by hereafter. The optimized station-specific measurements and their uncertainties given by Equations (5)-(8) are used as inputs for the inversion. To distinguish different events, we append a subscript to the notation of these inputs. For instance, has the same meaning as in Equation (5) but is for the th event. , , and are defined in a similar manner. Furthermore, the following vector notation is used: This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International The posterior probability of the model parameters ( ) can be written as where ( ) is the prior probability; ℒ ( , ) and ℒ ( , ) are the likelihoods of the time and amplitude data, respectively; and is a normalization constant. Direct computation of Equation (17) is infeasible because the normalization constant involves integration over the entire model space. However, the posterior probability can be estimated via probabilistic sampling, such as with the MCMC algorithm.
We assumed a Gaussian distribution for the prior probability of the horizontal locations, station correction terms, S-wave velocity, and attenuation factor: where and are the mean and standard deviation, respectively, and is either , , , , , or . We adopted Rayleigh distribution for event depths: Here, 0 is added to the usual formulation of the Rayleigh distribution. Without this term, the Rayleigh distribution is defined for positive values (i.e., > 0). Adding 0 changes the domain to > 0 . Introducing 0 may be useful for prohibiting hypocenters located above the seafloor, although we found that it did not affect the results significantly.  The likelihood function for the arrival time can be defined as follows: where is the predicted arrival time based on the hypocenter and S-wave velocity, and the subscripts and correspond to station and event indices, respectively. Because only the relative arrival times are known, Equation (20) includes an unknown event-specific time shift term . Rather than solving for this term, we set it to the MLE: Note that Equation (21)

corresponds to the averaged residual over stations weighted by
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International data variance. Similar equations hold for the logarithmic amplitudes: Notably, the term for source amplitude 0 is canceled out when Equation (12) is substituted into Equations (22) and (23), which eliminates the need to estimate the source amplitude beforehand.
Based on Equations (17)- (23), we can use the MCMC algorithm to sample the posterior probability. At each iteration, one of the model parameters is perturbed randomly, where the amount of perturbation is drawn from a zero-mean Gaussian distribution with a standard deviation, as given in Table 2. This modification is accepted or rejected in accordance with the Metropolis-Hasting criteria. In this study, we performed 4 million iterations, with the first 2 million iterations treated as a burn-in period. The sampled model parameters were saved at every 500 iterations during the second 2 million iterations.
We ran 100 chains of the MCMC algorithm in parallel and allowed them to mutually interact by using a parallel tempering technique for an efficient global search (Geyer, 1991;Sambridge, 2014). We set 20 Markov chains with a unit temperature, while the remaining 80 chains were assigned temperatures between 1 and 2000. A chain pair was randomly selected, and their temperatures were swapped at a certain probability to maintain the detailed balance (Sambridge, 2014). We repeated this temperature swap process 10 times per iteration.

Inversion results
We applied the above inversion method of Step 3 to the amplitude and time data from the 580 events that passed the quality control in Step 2. The likelihood almost monotonically increased with the number of iterations and converged within the burn-in period (Fig. 5a, black dots), which suggests that model parameters sampled after the burnin period can simulate the posterior probability. To evaluate the effect of the parallel tempering scheme, we conducted a parallel inversion analysis using 100 MCMC chains but without tempering. Although the likelihood increased at a slower pace than the tempered analysis (Fig. 5a, gray dots), most chains reached the same likelihood level as the tempering method at the 250,000th iteration (Fig. 5b). This good performance even without a tempering scheme likely reflects a well-posed inverse problem. As noted earlier, the quality control in Step 2 ensures a well-defined global minimum, which allows the inverse problem to be solved efficiently. This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International black and gray histograms show the results with and without the tempering scheme, respectively.
The inversion results are summarized in Fig. 6. Note that similar results were obtained from the non-tempered analysis (Fig. S2). The epicenters, which we defined as the mean of the MCMC samples, are tightly clustered in the map view. The typical horizontal location uncertainties (1 ) derived from the MCMC samples are <2 km in the east-west direction (blue histogram in Fig. 7a) and <3 km in the north-south direction (blue histogram in Fig. 7b). The errors are slightly less in the east-west direction than in the north-south direction because the seismic network geometry is elongated in the eastwest direction and variation of the subsurface structures is relatively gentle in the trenchparallel direction. The typical vertical uncertainties are < 5 km (blue histogram in Fig. 7c).
Unfortunately, the vertical uncertainties are insufficient to discuss the source faults of the tectonic tremors considering the subduction depth of ~6-8 km. Because of this loose constraint on the depth, some hypocenters are located above the seafloor. We may explicitly prohibit such unlikely solutions by increasing 0 in Equation (19), although this change had almost no influence on the horizontal locations (Fig. S3).
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International  The mean values of the correction factors range from -7.3 to 8.4 s for time delays and from -8.1 to 6.2 dB for amplification (Figs 6c and d). Overall, these values exhibit a smooth lateral variation, with stations near the trench experiencing earlier arrivals and a more significant amplification than predicted. The thinner accretionary prism near the trench likely explains the early arrivals, which allows seismic waves to travel through the subducted crust at faster velocities. In addition, the significant amplification at the trench is reasonable because the trench-fill sediments are less consolidated than the landward accretionary prism (Tsuji et al., 2011). Station MRE20 exceptionally shows a delayed arrival near the trench. Because this station is separated from the majority of events, this delayed arrival may account for the structural heterogeneities in the trench-parallel direction.
The posterior probabilities of the S-wave velocity and quality factor have narrow peaks, with mean values of 2.85 km/s and 274, respectively (Figs 6e and f). These values correspond to an attenuation strength of 2.00 × 10 -2 km -1 , and they are consistent with those obtained from the regression analysis (Fig. 4). The S-wave velocity of 2.85 km/s is somewhat slower than that reported for the oceanic crust of this region (>3 km/s) but is comparable to the velocity of the underthrust sediment immediately above the crust This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International (Akuhara et al., 2020). Yabe et al. (2021) independently estimated the attenuation strength of this region as a function of the hypocentral distance by using the seismic amplitudes of tectonic tremors that occurred in different periods, and their results are mostly consistent with our estimations.

Contributions of each factor
The proposed method offers several improvements compared to conventional analyses. For better understanding of its advantages, the contributions of different factors need to be considered, and hence we performed inversion under different settings (Fig. 8). Fig. 8(a) shows the inversion results from Fig. 6 (i.e., complete case). Fig. 8(b) shows the inversion results when the delay and amplification factors are excluded by setting their values to zero (i.e., without-correction case). Fig. 8(c) shows the inversion results when the time difference data are excluded and only the amplitude information was used (i.e., amplitude-only case). In this case, the amplification and quality factors are solved while the S-wave velocity is fixed at 3.0 km/s. This fixed value affects the estimation of through Equation (11) but not the other parameters. Fig. 8(d) shows the inversion results when only the time data are used with the S-wave velocity and delay factors solved (i.e., time-only case). This case does not involve the attenuation parameter.
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International In the without-correction case, the hypocenters are systemically located further seaward than in the complete case. Although we do not know the true hypocenters, the without-correction case shifts many events seaward of the trench, which is highly unlikely.
We conjecture that adding correction factors accounts for structural heterogeneities in the along-dip direction, which helps correct this artificial shift. The seaward shift is ~10 km on the western side, where station coverage is relatively limited. These shifts are greater than the uncertainties of the hypocenters shown in Figs 7(a) and (b). Failing to consider This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International these corrections can significantly bias the results and lead to misinterpretation.
The time-only case suffers from a greater uncertainty for the hypocenters than the amplitude-only case (Fig. 7). The hypocenters are more scattered in space in the time-only case (Fig. 8d) whereas they are similar to the complete case in the amplitude-only case (Fig. 8c). These discrepancies can be attributed to the considerable uncertainty in the time data, which can be quantitatively understood from the distance plots in Figs 3(c) and 3(f).
For example, the typical error in time can be read as 5 s from Fig. 3(c) while the typical error in amplitude can be read as 0.05 from Fig. 3(f). The typical errors in epicenter can then be calculated as and ⁄ represent the slopes of the regression lines in the time-distance and amplitude-distance plots, respectively. If ⁄ is 0.3 s ⋅ km −1 and ⁄ is 0.03, the error for the epicenter is 17 km using time data and 1.7 km using the amplitude data, which indicates a difference of an order of magnitude.
The large uncertainties in the time data can be attributed to the inconsistencies in the lag-time measurements between station pairs. Takemura et al. (2020) showed that a slow and heterogeneous accretionary prism complicates tremor waveforms as they propagate over longer distances. Lag-time measurements between stations at greater distances from each other are more susceptible to this waveform distortion, which can increase measurement inconsistency. A common strategy to mitigate this issue is to limit station pairs to those with shorter distances or high coherencies. However, such data selection is often based on subjective criteria.
Our results demonstrate the superiority of amplitude data for tectonic tremor location because it can pose tight constraints on hypocenters without any ad hoc selection of data. Challenges associated with using amplitude data may include difficulties with estimating the source amplitude, attenuation structure, and local site effects beforehand.
However, the proposed inversion approach eliminates the need for these prerequisite processes.

Spatiotemporal evolution of tremors
The proposed method provides well-constrained epicenters with typical estimation uncertainties of <2 km. This allows the spatiotemporal evolution of the tremor activity to be discussed in detail. As shown in Fig. 9(a), the located tremors can be divided into three main groups (A-C) separated by gaps of ~10 km. Fig. 9(b) shows that, during the first 9 days of the study period, source locations within group A mainly migrate to the east along the trough margin at a speed of ~2 km/day (phases i-iii). Such migration of tremors has been commonly observed worldwide, and it is thought to reflect an undergoing slow slip This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International event. In the middle of this period (i.e., phase ii), group B showed some tremor activity to the west. Although the gap between groups A and B makes it difficult to conclude the migration pattern, the tremor epicenters appear to show a westward migration at a speed of 8 km/s in phases i-ii. Once the main eastward migration reached group C, the migration speed abruptly changed to ~13 km/day (phase iv). This fast migration lasted for 3 days.
Then, after a quiescence period of about 1 day, relatively small-scale activity occurs in the northeastern part of Group A (phase v; see orange epicenters in Fig. 9a). The observed spatiotemporal evolution of the tremors is roughly consistent with that described by Ogiso & Tamaribuchi (2022), who used amplitude data from DONET stations to determine tremor locations. Temporal evolution of tremors projected along the X-Y profile (red line in (a)). The color notation corresponds to that in (a). (c) Probability of at least one tremor epicenter being This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International located within a 1 km × 1 km cell. The green dashed line represents the trench.
Our use of the stochastic sampling technique facilitated the exploration of subtle features within the tremor patterns while minimizing the risk of misinterpretation. For example, we can calculate the probability that any tremor epicenter is located at a particular geographical point ( , ) as , where ( , ) denotes the marginalized posterior probability for the th event epicenter ( , ). Visualizing this probability allows us to identify fine-scale spatial patterns of tremors without being disturbed by events with large uncertainties because such uncertain events have a limited impact on . Japan (Ide, 2010). The trench-parallel features of group C are likely associated with the topography of the decollement (Hashimoto et al., 2022). This disparity in the distribution patterns may correspond to differences in the migration speed, although a more detailed analysis is outside the scope of this study.

Conclusions and future perspectives
We proposed a novel three-step method for locating tectonic tremors that employs the optimization of pair-difference data, quality control via rough estimates of the propagation speed and attenuation strength, and joint inversion of multiple events using the MCMC algorithm. The proposed method eliminates the need for subjective tuning of data weights and avoids relying on prior knowledge of subsurface structures, local site effects, and source amplitudes. Although some subjective choices are still necessary to set quality control thresholds for and , these choices do not distort the uncertainty estimation. When applied to real data, the proposed method demonstrated its effectiveness.
Appropriately weighting data by their uncertainties was shown to mitigate the undesirable influence of low-quality data (Figs 8c and d), and the correction terms significantly reduced systematic biases (Fig. 8b). Furthermore, using a probabilistic mapping technique allowed us to better comprehend the detailed patterns in locations of tectonic tremors (Fig.   9c). Specifically, we were able to identify striations in the tremor sources. This provides valuable insights into the underlying structural factors that favor tremor activities.
This manuscript is a non-peer reviewed preprint submitted to EarthArXiv and is under consideration at Geophysical Journal International The proposed method still has room for improvement. One of the main assumptions is that the subsurface structures for and are uniform, which can potentially impact the results. Using more sophisticated correction factors, such as source-specific corrections, may help address this bias (e.g., Lomax & Savvaidis, 2022;Richards-Dinger & Shearer, 2000). Alternatively, the spatial variation of and can be solved as unknown parameters, similar to a tomographic approach. The narrow peaks observed in the posterior probabilities (Figs 6e and f) suggest that such an attempt could be promising.
One aspect that we did not discuss in the present study is the criteria for detecting tectonic tremors. In this study, we used the 98th percentile of the histograms of crosscorrelation coefficients as a threshold, which was an arbitrary choice. However, the propagation-based quality control in Step 2 of the proposed method provides an alternative approach to detecting tremors. Specifically, applying the selection criteria based on and to all time windows not prescreened by cross-correlation coefficients can incorporate wave-propagation information into the detection process, which would increase its robustness compared to relying solely on waveform similarities.
However, one drawback of this wave-propagation-based detection is that it requires high signal-to-noise ratios across the entire seismic network. Solving this problem is left for future work, but using such an objective detection method would help illuminate other important aspects of tectonic tremors, such as the frequency distribution (e.g., Nakano et al., 2019).
While obtaining the detailed features of tremor locations is key to understanding the physical processes behind them, it is particularly challenging for offshore regions, where the accurate location of tremors is hindered by strong heterogeneities in the shallow sedimentary structure. Our results demonstrated that our proposed method is applicable even to such challenging ocean-bottom data. Tectonic tremors that occur in shallow subduction zones remain underexplored. We believe that applying our proposed technique can shed new light on these phenomena.