Imaging different components of a tectonic tremor sequence in southwestern Japan using an automatic statistical detection and location method

In this study, we demonstrate the capability of an automatic network-based detection and location method to extract and analyse different components of tectonic tremor activity by analysing a 9-day energetic tectonic tremor sequence occurring at the downdip extension of the subducting slab in southwestern Japan. The applied method exploits the coherency of multiscale, frequency-selective characteristics of non-stationary signals recorded across the seismic network. Use of different characteristic functions, in the signal processing step of the method, allows to extract and locate the sources of short-duration impulsive signal transients associated with low-frequency earthquakes and of longer-duration energy transients during the tectonic tremor sequence. Frequency-dependent characteristic functions, based on higher-order statistics’ properties of the seismic signals, are used for the detection and location of low-frequency earthquakes. This allows extracting a more complete ( ∼ 6.5 times more events) and time-resolved catalogue of low-frequency earthquakes than the routine catalogue provided by the Japan Meteorological Agency. As such, this catalogue allows resolving the space–time evolution of the low-frequency earthquakes activity in great detail, unravelling spatial and temporal clustering, modulation in response to tide, and different


I N T RO D U C T I O N
Slow earthquakes refer to a class of events, corresponding to internal slow energy-release transients, that generate seismic or geodetic non-stationary signals with different characteristic timescales (e.g.Ide et al. 2007a;Obara & Kato 2016).Geodetic slow earthquakes (Hirose et al. 1999;Dragert et al. 2001;Lowry et al. 2001) are either long-term, with durations of months to years, or short-term, with durations of days, slow-slip transients.Seismic slow earthquakes span a wide range of characteristic durations and amplitude scales (Ide et al. 2007a;Peng & Gomberg 2010), producing non-stationary signals that are relatively depleted in high frequencies and different from those of 'regular' earthquakes.Weak continuous vibrations with no clear impulsive phases, known as tectonic tremor (Obara 2002), are among the most commonly observed seismic slow earthquakes signals.Tectonic tremor episodes typically last for several days and are often associated with aseismic short-term slow slip events, for example, 'episodic tremor and slip' (ETS; Rogers & Dragert 2003).Other seismic signals belonging to the slow earthquakes family and characterized by source durations of less than one second and few tens of seconds have been classified as low-frequency earthquakes (Katsumata & Kamaya 2003) and very-low frequency earthquakes (Ito et al. 2007;Ghosh et al. 2015), respectively.
Tectonic tremors (Fig. 1a) were first discovered and described more than 10 yr ago in the Nankai subduction zone, southwestern Japan, by Obara (2002).Since then, both geodetic and seismic slow earthquakes have been observed in a wide variety of subduction settings (Schwartz & Rokosky 2007;Peng & Gomberg 2010;Obara 2011;Ide 2012) and in the vicinity of intracontinental fault zones (Nadeau & Dolenc 2005).Analysis of well recorded tectonic tremor episodes in several regions pointed out that a significant fraction of their signals appears as mixing of impulsive low-frequency earthquakes (Shelly et al. 2007a;Brown et al. 2009).However, the mechanisms of this mixing are not yet fully understood.Moreover, spectral similarity between the two types of signals and consistency of low-frequency earthquakes focal mechanism with shear slip along the plate interface (e.g.Shelly et al. 2007a;Ide et al. 2007b) led to the current view of tectonic tremor being a swarm of low-frequency earthquakes.While the physics of their source mechanism is not fully understood, slow earthquakes (both seismic and geodetic) are often considered to be related to fluid-enabled shear transients at brittle-ductile transition zones of active faults, and associated with slow transient energy release processes that characterize slow dynamics of faults driven far from equilibrium, during which the build-up of mechanical energy may sometimes lead to a major earthquake (e.g.Heinrich et al. 2010;Benzi et al. 2016;Obara & Kato 2016).As such, slow earthquakes potentially represent a unique source of information for understanding the physics of the 'preparation' processes of these large events and for improving monitoring and forecasting methods (Brudzinski 2011;Obara & Kato 2016).
Detection, location and source restoration of tectonic tremor episodes challenge traditional methods used for the analysis of 'regular' earthquakes.In contrast to the latter, seismic signals associated with tectonic tremor are relatively weak, low-amplitude and longer-duration non-stationary signals (with no definite beginning), sometimes similar to ambient noise.Their main characteristic feature (e.g.Obara 2002) is a long-duration transient energy increase recorded simultaneously at multiple stations (Fig. 1b).
Since the work of Obara (2002), todays' standard techniques for detection and location of tectonic tremor, are based on 'network response function' (NRF) such as amplitude-based NRF (Nadeau & Guilhem 2009;Husker et al. 2010Husker et al. , 2012) ) or beamforming-like techniques (e.g.back-projection, reverse migration, source scanning) exploiting traveltime differences between the network receivers estimated from cross-correlations of smoothed, averaged energy envelopes in a grid-search procedure, and, in some cases, combined with energy attenuation analysis (Obara 2002;Kao & Shan 2004;Kao et al. 2005;La Rocca et al. 2008;Payero et al. 2008;Wech & Craeger 2008;Ghosh et al. 2009Ghosh et al. , 2010a;;Maeda & Obara 2009;Suda et al. 2009;Ide 2010Ide , 2012;;Sun et al. 2015).Energy envelope functions used in these analyses are typically calculated from single or three component seismograms band-pass filtered in the 2-10 Hz frequency band, corresponding to the observed predominant frequency range of tremor signals.For most of these techniques, resolution of the tectonic tremor source locations remains limited.Even though, careful strategy for detection and removal of outliers, together with spatial clustering analysis (Wech & Creager 2008;Ide 2010), may reduce errors of tremor source location to less than 5 km, these are post-processing schemes that are not data-driven.
Very few studies on detection, location and characterization of tectonic tremor attempted to exploit directly the phase information of the recorded waveforms.The recent work of Cruz-Atienza et al. (2015) provided a high-resolution scheme for simultaneous location and characterization of tectonic tremor sources in Guerrero, Mexico, by combining the amplitude method of Husker et al. (2012) and a waveform polarization analysis.Although very promising, the method requires preliminary identification of the periods of tectonic tremor activity either by visual inspection of the records or specific automated scheme (e.g.Brudzinski & Allen 2007;Husker et al. 2010;Sit et al. 2012).
There currently is a number of available large tectonic tremor catalogues the catalogue of the National Institute for Earth Science and Disaster Prevention (NIED), providing hourly location of tectonic tremor sources in southwestern Japan using the method of Maeda & Obara (2009) and the clustering scheme of Obara et al. (2010); the more recent catalogue by Annoura et al. (2016) for southwest Japan, using refined procedure for the extraction of tectonic tremor sequences (periods of continuous energetic tectonic tremor activity); the Automatic Tremor Monitoring System catalogue of Suda et al. (2009); the catalogue of tectonic tremor activity in Cascadia, using the method of Wech & Creager (2008); and the world tectonic tremor database (Idehara et al. 2014), based on the method of Ide (2010Ide ( , 2012)).There also currently exists an ongoing effort in building a centralized database combining different catalogues in a unified format in a single repository at the University of Tokyo (Kano et al. 2017).
A different approach of analysing tectonic tremor activity, after identification of the periods of tremor activity, focuses on the short-duration impulsive transients, such as low-frequency earthquakes, that are mixed within the tectonic tremor signals.For some of those, the S-phase (and sometimes P-phase) arrivals can be manually picked and the corresponding sources can then be located using traditional earthquake location methods.The Japan Meteorological Agency (JMA) includes such manually located low-frequency earthquakes into its routine catalogue, distinguishing them from 'regular' earthquakes (Katsumata & Kamaya 2003).The extracted waveforms of individual low-frequency earthquakes can be later used to derive and statistically analyse their source characteristics (Shelly et al. 2006;Bostock et al. 2015Bostock et al. , 2017;;Thomas et al. 2016).
Low-frequency earthquakes can also be detected and extracted from low signal-to-noise ratio (SNR) time-series using matchedfilter techniques, exploiting their impulsive signature and repetitive behaviour (e.g.Shelly et al. 2006Shelly et al. , 2007a;;Shelly & Hardebeck 2010;Bostock et al. 2012;Frank et al. 2013;Frank & Shapiro 2014).Such template-based methods provide more complete catalogues, with high temporal resolution, but require a rather comprehensive initial dictionary of template events and significant computational resources.Although powerful, template-based methods are intrinsically limited by the assumption that low-frequency earthquakes activity can be fully represented by a set of repeating events.Up to now, alternative data-driven detection and location methods are very few: for example, subspace detector method by Maceira et al. (2010) and intercomponent cross-correlation of La Rocca et al. (2009).Both studies, however, are restricted to the analysis of limited data sets of 1 hr and 2 weeks respectively.Their applicability and performance for more general case remains to be evaluated.
The nature of the mixing interactions of low-frequency earthquakes, and the relation between these impulsive short-duration transients and longer-duration energy transient components during tectonic tremor sequences is still not fully understood.Accurate detection and source location of these two signal components is important for restoring the energy-radiation source regions and gaining new insights into low-frequency earthquake activity and multiscale interactions during the tectonic tremor sequences.Therefore, there exists a need in developing and evaluating, efficient automatic detection and location methods free from a priory assumptions or parametrization.
In this study, we demonstrate the capability of the recently developed network-based BackTrackBB method (Poiata et al. 2016) to analyse short-duration impulsive transients associated with lowfrequency earthquakes and longer-duration energy transient components during a 9-day long (from 2012 May 25 to June 2) energetic tectonic tremor sequence in southwestern Japan.The paper is organized as follow: Section 2 describes the study area and the analysed data.Section 3 provides a brief overview of the BackTrackBB method.In Section 4, the BackTrackBB methodological details relevant to the detection and location of low-frequency earthquakes during the tectonic tremor sequence are presented.The performance of the method is illustrated by comparing the derived catalogue of low-frequency earthquakes with the JMA catalogue for the same tectonic tremor sequence.Then, the space-time activity of the lowfrequency earthquakes is analysed in more detail.In Section 5, the BackTrackBB methodological details relevant to the detection and source location of longer-duration energy transient components during the same tectonic tremor sequence are presented.We derive a catalogue of longer-duration energy sources characterized by their duration and 3-D spatial likelihood maps of the emitting source areas.The summary 3-D likelihood map of the 9-day tectonic tremor sequence, built from this catalogue, is then analysed with respect to the extracted low-frequency earthquakes activity.
We then examine the two catalogues analysing the spatial relation between the extracted short-duration low-frequency earthquakes activity and the longer-duration energy sources.This analysis provides a potential constrain on the size of the energy-release source regions in relation with low-frequency earthquakes clustering activity during the analysed tectonic tremor sequence.
We analysed a tectonic tremor sequence that occurred beneath western Shikoku in southwestern Japan (Fig. 1a).This area is known as having one of the highest rates of tectonic tremor seismic energy release along the Nankai subduction zone (e.g.Yabe & Ide 2014;Annoura et al. 2016).Very-low-frequency earthquakes (Ito et al. 2007) and slow-slip events, have been also observed in this area (Hirose & Obara 2005).ETS in Nankai subduction zone are located in the transition part of the subduction slab, both shallower and deeper than the seismogenic areas of the past megathrust earthquakes (Yamashita et al. 2015;Obara & Kato 2016).While deep tectonic tremors and low-frequency earthquakes, have been extensively studied since their discovery (Obara 2002), it is only recently, that shallow tectonic tremor was observed and analysed (Yamashita et al. 2015).Sources of deep tectonic tremor are known to be continuously distributed within a narrow belt-like zone along the downdip edge of the megathrust seismogenic region at depth of around 35 km (e.g.Obara & Kato 2016).Analysis of the impulsive short-duration transients, such as low-frequency earthquakes, together with longer-duration energy transient components during a sequence of such deep tectonic tremor activity represents the main goal of our study.
We used velocity seismograms recorded at 21 stations of the highsensitivity borehole seismic network (Hi-net) in western Shikoku, maintained by NIED, Japan (Okada et al. 2004;Obara et al. 2005; Fig. 1a).The continuous 9-day long seismograms cover the time period of 2012 May 25-June 2 Japan Standard Time (JST), corresponding to an energetic tectonic tremor sequence.Although Hi-net stations are equipped with three component sensors, owing to the predominance of the S-wave component in recorded tectonic tremor signals pointed out by numerous studies (e.g.Obara 2002;Shelly et al. 2006;Maeda & Obara 2009), only the NS and EW horizontal components have been used here.
JMA differentiates since 1999 low-frequency earthquakes (Fig. 1c) as a separate class of events in its seismic catalogue (Katsumata & Kamaya 2003).Low-frequency earthquakes reported in this catalogue are located, following the standard location scheme routinely applied for 'regular' earthquakes, using manually picked S-and, when possible, P-wave.In this study, we used the JMA catalogue of low-frequency earthquakes for validating our results.

G E N E R A L OV E RV I E W O F T H E M E T H O D
BackTrackBB (Poiata et al. 2016) is an automatic detection and location method that exploits the coherency of multiscale, frequencydependent characteristics of non-stationary seismic signals recorded across a network of stations.The workflow of the method can be, in general lines, separated in two phases: signal processing and detection and location.The signal processing part builds time-series of characteristic functions (CFs) by extracting targeted properties of the continuous non-stationary signals.This step takes into account the potential multiscale nature of the signal's characteristics and incorporates time-frequency analysis to extract relevant frequency-dependent information.The time-series of CFs are then used in the detection and location process, which, itself, consists of two parts: the station-pair time-delay functions estimation, exploiting information on the coherency of the CFs between pairs of stations, and the spatial mapping and stacking of station-pair timedelay estimate functions.This leads to time-series of 3-D spatial likelihood images of source locations on which detection-location is performed.Each phase of the method can be specialized and tailored depending on targeted features of the signal, such as impulsive short-duration transients, like low-frequency earthquake signals, or longer-duration energy transient components during tectonic tremor.This distinguishes BackTrackBB from other, similar, existing methods, such as Source Scanning Algorithm (Kao & Shan 2004), Waveloc (Langet et al. 2014) or Loki (Grigoli et al. 2014).In particular, with time-localized CFs based on frequency-dependent, higher-order statistics (HOS) representation of the recorded signals, BackTrackBB can automatically and efficiently detect and locate impulsive, short-duration low-frequency earthquakes during tectonic tremor activity.In contrast to matched-filter techniques, the method does not require building a dictionary of template events while providing a more complete and better time-resolved event catalogue than the catalogue of JMA.This high-resolution catalogue allows analysing in great detail the space-time activity of low-frequency earthquakes during tectonic tremor sequence and resolving migration patterns, spatial clustering and temporal modulation.On the other hand, with CFs based on smoothed frequencydependent energy envelopes, BackTrackBB can exploit coherency of longer-duration energy transients in the signal and provides information on extended energy-release source area during tectonic tremor sequence, in contrast to classical time and space-averaged point source location (Maeda & Obara 2009;Suda et al. 2009;Ide 2010).

Methodology
In this section, some aspects related to the different steps of the BackTrackBB method (Poiata et al. 2016) are discussed in the context of the automatic detection and location of short-duration, impulsive low-frequency earthquakes during a tectonic tremor sequence.

Signal processing scheme: multiband statistical signal representation
The signal processing phase transforms the non-stationary signal recorded at each station into a characteristic function (CF) designed to extract and localize in time the time-frequency short-duration (impulsive) transient features that are of most relevance for the signature of low-frequency earthquakes (Figs 1b and c).For such short-duration impulsive signals, the CF is constructed using higherorder signal statistics (HOS-CF) that efficiently capture and extracts changes in statistical properties-not biased by relative amplitudeassociated with short impulsive transients in complex signals with low SNR (Saragiotis et al. 2002;Lokajíček & Klíma 2008;Langet et al. 2014).
The HOS-CF and, correspondingly, the HOS moments of the signal can be efficiently calculated using a recursive scheme based on exponential time averaging as where μ, is the first central moment (mean) expressed as μ ; m2 , mn are estimates of the second (variance) and nth central moments of the signal, respectively; u(t i ) denotes the ith-sample of the signal; and C: 0 ≤ 1 − C ≤ 1 is the decay constant parameter, C = T /T decay , associated with the characteristic time decay, T decay , of the exponential averaging with T being the data sampling interval.Although, HOS of any order n ≥ 4 can be used, in this study, we applied fourth-order statistics leading to the so-called kurtosis HOS-CF.
In order to take into account time-frequency variability of lowfrequency earthquakes signals, a frequency-dependent strategy is adopted for the construction of the broad-band HOS-CF (CF HOS ) based on multiband filter (MBF) decomposition of the signal (Lomax et al. 2012).According to the MBF algorithm, the original signal is run through a pre-defined bank of narrow band-pass filters, with n = 0,. . ., N band central frequencies, covering the frequency interval of interest [ f min , f max ] for low-frequency earthquakes signals.Frequency-dependent HOS-CFs, CF(t, f ), are constructed for each of the narrow-band filtered signals using eq.( 1).The final broad-band CF is then computed using the maximum operator as follows: (2) This time-frequency decomposition approach allows retaining at each time the strongest HOS signature of short-duration, impulsive transients in the different frequency bands.It also provides an efficient frequency-dependent denoising scheme in the case of low SNR.Finally, to improve the time-localization of the impulsive transients, the positive time derivative of the broad-band CF HOS is convolved with a Gaussian window of half-width .σ = T decay /2, providing an a priori uncertainty measure associated with HOS evaluation.As a result, the final CF signature of a transient is a symmetric Gaussian, centred on the onset of the transient (e.g.phase arrival; Figs 2e and 3a).
For the HOS-CF recursive evaluation scheme (eq.1), when focusing on the signals of short-duration, impulsive low-frequency earthquakes, a constant and frequency-independent T decay is selected and set to 2.0 s (guided by the predominant source duration of lowfrequency earthquakes ∼1 s).The applied MBF filter-bank is composed of 20 logarithmically spaced filters covering the frequency range of 2.0-15.0Hz, adapted to the predominant frequency range of low-frequency earthquakes.Processing of the 9-day long continuous signals, recorded at selected Hi-net stations, is implemented by segmenting the signal into 1 hr time-windows with an overlap of 60 s at the ends.An example of the MBF signal decomposition and frequency-dependent HOS-CF is provided in Figs 2(a)-(d), for a selected low-frequency earthquake extracted from the JMA catalogue and recorded at one of the Hi-net stations.The final broadband CF HOS is shown in Figs 2(d) and (e).It illustrates how the time-frequency decomposition helps to enhance the S-wave phases associated with the low-frequency earthquake and to reduce the contribution of accidental spikes, such as the one at around 1290 s in Fig. 2(e).

Time-delay estimation between pair of stations
The second phase of BackTrackBB performs the evaluation of the station-pair, time-delay estimate (TDE) distributions (functions) extracting the information about the coherency of CFs between the pairs of stations.Details regarding this step can be found in Poiata et al. (2016) and only relevant aspects for this study are discussed here.
The calculation of the TDE distributions is based on local crosscorrelation (LCC), implemented with the efficient recursive algo-rithm of Hale (2006).The LCC is defined in the time-domain and is computed at each time-sample within a localized time neighbourhood, providing an estimation of the time-dependent similarity between two signals.Conceptually, this is equivalent to constructing before the cross-correlation windowed functions, associated with each station-pair, and estimated from the final CFs f (t) and g(t), by convolving them with a smooth Gaussian window: where is the Gaussian windowing function of half-width σ , centered on t = t .
Using the standard formulation of cross-correlation, and denoting by l the cross-correlation time-lag, the LCC can be expressed after the change of variable where v(t, l) is a Gaussian windowing function defined now in the time-lag space as 4σ 2 e The 2-D LCC function defined in eq. ( 4) provides the TDE distribution at times t = t , while the parameter σ of the 2-D Gaussian filtering (eq.5) controls the scale of the temporal coherency retained in the LCC function and, correspondingly, in the TDE distribution.
Although, the TDE calculation step is conceptually independent from the choice of the CF, the parametrization needs to be adapted to the expected characteristics of seismic signal.For impulsive shortduration transients such as low-frequency earthquakes, σ must be chosen so, as to capture the coherent information between the CFs of a station-pair at short time-scales.The value of σ = 0.5 s has been chosen in this study.
It is worth noting that, as a result of the signal-processing step, the broad-band HOS-CF is Gaussian-shaped by construction and is centred at the onset of the transient (phase arrival; Fig. 3a).Therefore, increasing σ in the 2-D Gaussian filtering results in increasingly smoothing and broadening the TDE distributions, which in the following detection and location step will be transformed-through an imaging function-into a wider spatial likelihood of the radiatingsource volume.This effect of the Gaussian filtering is illustrated, for one selected low-frequency earthquake and different values of σ , in the Supporting Information (Fig. S1).
For low-frequency earthquakes analysis, the LCC-based estimation of the TDE distributions is performed on 1 hr segments of the 9-day continuous broad-band HOS-CFs, using sliding timewindows of 20 s (Fig. 3a) and an overlap of 18 s between two consecutive windows.The CFs are downsampled to 20 Hz which significantly reduces the LCC computation times.The final timedependent function of the TDE distribution, used for the detection and location phase, is then coarse-grained according to the size w of the sliding time-window as The size of the sliding window is selected taking into account the maximum interstation distance among the station pairs, set here to 100 km, and considering the short duration of low-frequency earthquakes.

Detection and location scheme
The detection and location phase of BackTrackBB is based on mapping (projecting) the station-pair TDE distribution for each time-window (eq.6), onto a 3-D spatial grid of theoretical traveltime differences calculated under the assumptions of a velocity model and a seismic phase (e.g.S-phase here).Summation of the projected TDE distributions over all the set of station-pairs leads to time-series of 3-D grids that provide the spatial (normalized) likelihood for the source to be located at a given position (Fig. 3b).
In this study, a 3-D grid covering the region of 292 km × 190 km × 100 km in western Shikoku (Fig. 3b), with a 1 km 3 element grid-size, has been defined and S-wave was assumed to be the predominant phase (on the horizontal component) of the tectonic tremor signals (e.g.Obara 2002;Shelly et al. 2006).Theoretical S-wave traveltimes, associated with the 3-D grid, were estimated using the Grid2Time routine of NonLinLoc (Lomax 2005(Lomax , 2008)), and the 1-D velocity model of Kubo et al. (2002).
The BackTrackBB detection and location step is carried out for each of the 20 s sliding time-windows (with an overlap of 18 s), used for the estimations of TDE distributions.The selection of the stations used for detection and location of low-frequency earthquakes is an important issue as they are characterized by (relatively) weak signals (low SNR), often mixed with background tectonic tremor signal.A heuristic assumption in this study is that low-frequency earthquakes are detectable at stations within 50 km epicentral distance from the source (e.g.Fig. 1c).Accordingly, the selected 21 stations in western Shikoku are divided into subgroups of 13-17 stations; each consisting of stations located within 50 km distance from the mean NIED daily location of tectonic tremor sources.In general, this station-selection criterion is well adapted to the expected slow migration pattern (∼10 km d −1 ) and the clustering of the low-frequency earthquakes previously reported in this area (e.g.Obara 2002;Shelly et al. 2006Shelly et al. , 2007a)).
Short-duration, impulsive signals of low-frequency earthquakes are assumed to be associated with a point source.Accordingly, detection of a low-frequency earthquake is declared in this study each time the normalized value of the final 3-D spatial likelihood map (imaging function) exceeds the 0.65 (65 per cent) threshold.This value has been selected after a preliminary calibration analysis, where a 1 hr record of tectonic tremor activity was analysed with the BackTrackBB method and the resulting event catalogue compared with the JMA catalogue for the same period.The threshold value was set to maximize the detection of the low-frequency earthquakes already reported in the JMA catalogue and to minimize the number of false detections.All the detected events in this selected 1 hr data set went also through a visual inspection and validation.Example of a typical imaging function for a single detected low-frequency earthquake is shown in Fig. 3(b).The imaging function is well focused around a likelihood maximum that exceeds the prescribed threshold, triggering the detection.The event location is then associated with the space-time location (X, Y, Z, T 0 ) of the likelihood maximum.As shown in Fig. 3(c), the location of the detected lowfrequency earthquake is in very good agreement (within ∼5 km of the horizontal distance off-set) with the JMA hypocentre (marked as white star in Fig. 3c) that falls within an area of the imaging function corresponding to a likelihood greater or equal to 0.9.
The entire data set of the 9-day recordings, once processed, leads to a preliminary catalogue containing all the detected and located events, together with the origin times and the arrival times of the S-phases used for location (Fig. 3a).This catalogue undergoes a post-processing step that groups together multiple detections arising from the significant overlap (90 per cent) between the consecutive sliding time-windows and the move-out of the phase-arrival for the selected stations.When at least two consecutive detections provide locations within 2.0 km distance from one another and differences in origin time of less than 1.5 s, these detections are associated with a single event and the event location is assigned to the largest value of the maximum likelihood among the different imaging functions.

Results
The final catalogue of low-frequency earthquakes obtained for this 9-day long tremor sequence consists of 2212 detected and located events (Fig. 4).Investigation of the spatial and temporal distribution of the detected low-frequency earthquakes during this period, shows an activity that started in the southwestern part of western Shikoku and extended with time towards the northeast, along the strike direction of the subducting slab (Figs 4a and b), covering an area of about 120 km in length and 50 km in width (Fig. 4a).The majority of events are constrained in depth between 30 to 40 km as shown in Figs 4(b) and (c).The relative scattering of the depth locations comes, on one side, from the use of S-phase information only for the detection and location, and, on the other side, can be attributed to the fact that the catalogue is the direct output of the automatic scheme with no strict revision (based on events' waveform evaluation).Nevertheless, a careful inspection of the along-strike depth profile (Fig. 4b) indicates that the detected low-frequency earthquakes follow the depth-profile of the subducting slab, with the highest density of events being concentrated around the depth limit of 30-40 km.It is worth noting here that the hypocentre distribution of the detected events reflects the slight deepening of the slab towards southwest confirmed by the average depth profile fitted to the along-strike 4 km binned hypocentre distribution (Fig. 4b).The final catalogue also well resolves the evolution and the complexity of low-frequency earthquakes activity in space and time.During the first 3 day, events are relatively uniformly distributed in the southwestern area, while in the following days the low-frequency earthquakes activity exhibits clear spatial clustering in the central and north-eastern parts of the region.The most active clusters of low-frequency earthquakes are observed during the 4th-to-7th days (84-168 hr, Fig. 4a) of the sequence, around the along-dip depth profiles 2-2 and 3-3 in Fig. 4(c) and just above the 'mantle wedge' (Figs 4a and c).Such clustering of the activity, was previously reported in the same regions, and analysed by Shelly et al. (2007a,b) for selected (6 day long) tremor episodes and Ide (2010) for 5 yr continuous recordings of Hi-net stations.According to these studies, the clear along-strike clustering is likely to reflect some structural or geometrical complexity of the plate boundary.Similar observations were also reported in Cascadia (e.g.Ghosh et al. 2010b) and Guerrero, Mexico subduction zones (e.g.Peng & Rubin 2017) and interpreted as resulting from 'frozen' heterogeneities in the slowslip source regions.

Comparison with JMA low-frequency earthquakes catalogue
The low-frequency earthquakes catalogue has been compared with the existing catalogue provided by JMA for the same time-period (Fig. 5).The JMA catalogue provides 341 low-frequency earthquakes during this 9-day tectonic tremor sequence (Fig. 5c).Although the epicentres of the events from JMA cover mostly the same area and have similar large-scale features, the 2212 low-frequency earthquakes detected and located in this study (Fig. 5a) provide a more detailed and resolved picture of the spatial and temporal activity of low-frequency earthquakes during the tectonic tremor sequence (Fig. 5).The depth resolution of the low-frequency earthquakes from the JMA catalogue is more accurate than the one provided by the automatic BackTrackBB catalogue (Figs 5a and c), as the JMA procedure involves manual detection and picking of the S-phases (and in some cases of the P-phases).Nevertheless, a more detailed investigation of the 290 events (85 per cent of JMA catalogue) that are common to both catalogues (Supporting Information Figs S2 and S3) indicate that the automatically located lowfrequency earthquakes in this study are in good agreement with the locations of JMA.A remarkable feature resolved in the automatic BackTrackBB catalogue is the 24 hr modulation of the low-frequency earthquakes activity likely associated with the tidal force, as shown by the hourly activity rate in Fig. 5(b).An additional inspection of long-duration waveforms (>20 h) recorded at multiple stations and spectrograms supports this observation, ruling out the influence of the day versus night cultural noise.The JMA catalogue does not resolve this feature of the tidal influence (Fig. 5d).Tidal effect is more evident around the main clustering activity period, 4th-to-7th days (84-168 hr), of the tectonic tremor sequence.Tidal triggering-and-modulation of tectonic tremor and low-frequency earthquakes is well documented by previous studies in this area of western Shikoku (Shelly et al. 2007b;Ide 2010) as well as in other areas, such as Cascadia (Rubinstein et al. 2008;Houston et al. 2011;Thomas et al. 2013).However, this is the first time that such an observation in western Shikoku can be resolved during a single and relatively short duration (9-day) tectonic tremor sequence thanks to the rate of low-frequency earthquakes detection provided by the direct and automatic analysis.
A close comparison between the JMA and the automatic Back-TrackBB low-frequency earthquake catalogues indicates that both capture similar features of activity during this tectonic tremor sequence: along strike propagation, spatial and temporal clustering.The space and time location of the main clustering activity of lowfrequency earthquakes extracted from the two catalogues is very similar.Such clusters of intense activity might correspond to source volumes that generate energetic low-frequency earthquakes, which can be detected and located manually in the JMA catalogue (see Shelly et al. 2007b).Nevertheless, the automatic BackTrackBB catalogue is significantly richer and more resolved in time, with ∼6.5 times more events providing details of the space-time lowfrequency earthquake activity during the 9-day tectonic tremor sequence.

Space-time evolution of low-frequency earthquakes
Three well-resolved space-time patterns of the low-frequency earthquakes activity could be extracted from the analysis of the automatic BackTrackBB low-frequency earthquake catalogue.
The first pattern is a large-scale, slow along-strike propagation of the main activity-front over distances of ∼120 km, as observed in Figs 4(a) and 5(a).Plotting the along-strike distribution of the detected low-frequency earthquakes as a function of time (Fig. 6a), an averaged propagation speed of the main activity-front is estimated to be 10-12 km d −1 .This is in good agreement with previous studies reporting such slow (5-20 km d −1 ) along-strike migration as one of the main features of the tectonic tremor activity in Shikoku, Japan (Obara 2002;Obara & Sekine 2009), in Cascadia (Dragert et al. 2004;Kao et al. 2006;Wech et al. 2009;Houston et al. 2011), and more recently in Guerrero, Mexico (Peng & Rubin 2017).
The two other patterns are associated with more complex short-scale migrations superimposed on the large-scale slow migration process, and having 10-100 times faster characteristic speeds.
The second space-time pattern corresponds to fast (5-7 km hr −1 ) short-scale along-strike migration occurring behind the main activity-front of the low-frequency earthquakes, and propagating into the backward direction.This pattern is observed as elongated stripes in Fig. 6(a), especially clear during the most active clusters of low-frequency earthquakes marked by red-shadowed rectangles (1-3) in Figs 6(a) and (b), that are also associated with high tectonic tremor signals' energy (Fig. 6c).Huston et al. (2011) first reported such rapid reverse-strike activity migrations-the so-called rapid tremor reversals-in Cascadia.This study provides the first evidence of such rapid low-frequency earthquakes activity reversals in western Shikoku.It is worth noting here that the clearest reverse-strike propagation patterns are observed during periods of intense low-frequency earthquakes activity (Fig. 6b) that are also associated with high tectonic tremor signals' energy (Fig. 6c).They initiate at the large-scale migration front and propagate in opposite direction-towards southwest-over distances of about 20-25 km.Four clear episodes of reverse-strike propagation are marked by arrows in Fig. 6(a), which compares the along-strike time distribution of low-frequency earthquakes activity with the hourly rate of detected events and the energy envelopes of seismograms at three selected stations (Fig. 6c).A zoom on the strongest cluster (number 3) provides more details about the propagation pattern (Fig. 7), allowing the identification of the migration onset and the estimation of the propagation speed at ∼5 km hr −1 (Fig. 7b).This episode is actually associated with high amplitudes of the tectonic tremor signals as shown in Fig. 8.These observations are in good agreement with previous results in Cascadia (Houston et al. 2011;Thomas et al. 2013), even though, the reported speed of 5-7 km hr −1 in this study is slightly slower than the speeds of 7-17 km hr −1 estimated there (e.g.Houston et al. 2011).Whether this reverse-migration speed does represent a characteristic feature of the Shikoku subduction zone is currently unclear, and further studies, systematically analysing the tectonic tremor sequences in the Shikoku area, are needed to answer this question.
The third (last) resolved space-time pattern corresponds to rapid along-dip migrations.It can be observed during the most active clusters of low-frequency earthquakes as 15-20 km long, subparallel stripes superimposed on the large-scale pattern of low-frequency earthquakes activity (Fig. 6b).A zoom on one of these clusters (Fig. 7) shows the complexity of these short-scale along-dip patterns, consisting of successive downdip and updip migrations with propagation speeds varying from 7-10 to 75-80 km hr −1 .This is up to 10 times faster than the observed reverse-strike migration speeds.In contrast to the latter, fast short-scale along-dip migrations do not seem to be related to the most active period of the clusters nor to the high amplitude of the tectonic tremor signals, and they can be observed through the entire cluster period.Such fast alongdip migrations were previously reported and studied in this area by Shelly et al. (2007a,b), while similar observations in Cascadia were reported from dense array beamforming analysis by Ghosh et al. (2010b) who observed 'streaks' of tremor propagating back and forth parallel to the plate convergence direction with speeds of 25-100 km hr −1 .More recently, Peng & Rubin (2017) reported similar features in Guerrero, Mexico.
The more complete and time-resolved catalogue of lowfrequency earthquakes provided by the BackTrackBB method during this single 9-day tectonic tremor sequence, allows extracting rich information about the space-time activity patterns in the western Shikoku area.These observations provide important information for improving our understanding of the driving physical processes controlling short-scale low-frequency earthquakes activity, possibly in relation with slow-slip and diffusion processes at larger space and timescales.However, such investigations are beyond the scope of this study.

S PA C E -T I M E I M A G I N G O F L O N G E R -D U R AT I O N E N E RG Y S O U RC E S D U R I N G T E C T O N I C T R E M O R
Seismic signals associated with tectonic tremor are relatively weak and low SNR long-duration transients with no well-defined beginning (Fig. 1b).The most easily observed feature of tremor signals is a long-duration transient energy increase recorded simultaneously at multiple stations (e.g.Obara 2002 and Fig. 1b).Even though low-frequency earthquakes are very valuable source of information, they only account for short-scale energy-release processes during tectonic tremor sequence.Detecting and locating sources of longerduration energy transient components is important for understanding their relationship with impulsive low-frequency earthquakes activity and multiscale interactions during the tectonic tremor sequences.
In the following, we present a first study based on the extension and specialization of the BackTrackBB method for the detection and location of longer-duration energy transients during a tectonic tremor sequence.The objective is to focus on longer-duration signal component that might be associated with possible low-frequency earthquake interactions during the tectonic tremor energy-release process.This is achieved by restoring longer-duration radiating sources as time-series of 3-D images of source likelihood location and analysing them in relation with the low-frequency earthquakes activity previously extracted.

Signal processing scheme
The characteristic functions, targeting detection of longer-duration energy transients during tectonic tremor sequence, are constructed from frequency-dependent energy envelopes.More precisely, we estimate the root-mean-square (RMS) energy envelope recursively as where C = T /T decay , with 0 ≤ 1 − C ≤ 1, is a constant associated with the characteristic time-decay, T decay , of the exponential averaging, and T is the data sampling interval.
The CF estimation is again based on MBF time-frequency signal decomposition accounting for possible frequency-dependent energy variations during the tectonic tremor sequence.The MBF signal decomposition is performed, for consistency, with the same filter-bank as in the HOS-based low-frequency earthquakes analysis.The time-frequency energy envelopes, CF env (t, f ), are calculated for each of the narrow-band filtered signals using eq.( 7), and the final broad-band CF is constructed using a summation operator as In contrast to the HOS-based CF, the decay constant T decay is assumed to be frequency-dependent.A value T decay = 0.75 s is chosen for the central frequency f N band/2 of the MBF filter-bank, while for other frequencies, T decay is estimated using the following scheme: where f n is the frequency of the nth narrow-band filter.The signalprocessing steps and the finally reconstructed broad-band CF are illustrated Fig. 9, for a 1 hr signal recorded at a single station.In practice, for the 21 Hi-net stations, the broad-band CFs are estimated from 1 hr segments of the 9-day continuous records, with an overlap of 60 s at the ends of each segment.

Time-delay estimation and detection and imaging scheme-reconstructing longer-duration energy-release sources during tectonic tremor sequence
Detection and location (imaging) of long-duration energy release sources follows similar steps as those for the analysis of the low-frequency earthquake previously discussed in Sections 4.1.2and 4.1.3,with some specific modifications that are presented below.
The analysis is performed over 1 hr segments of the energy envelope CFs using sliding time-windows of 180 s with an overlap of 140 s (∼78 per cent) between two consecutive windows.The CFs are down-sampled to 1 Hz sampling rate, in agreement with the longer-duration of the targeted signals (expected to be 2 s).The size of the sliding window is selected to account for the targeted longer-duration energy-release transients (of the order of minutes) within the non-stationary tectonic tremor signal.Previous studies related to tectonic tremor source location, based on energy-envelope cross-correlation techniques, typically use timewindow sizes from 60 s (Obara 2002;Maeda & Obara 2009) to 300 s (Ide 2010).
The detection and imaging step, first involving the LCC-base estimation of the station-pair TDE distributions (Section 4.1.3),is tuned now for extracting longer-duration coherency information from the CFs.Then, at each time, the normalized station-pair TDE distributions are projected onto a 3-D spatial grid using theoretical S-wave time-delay estimates, leading to likelihood images, that is, imaging functions, which are stacked over all station-pairs.The time-series of 3-D images provide information on the likelihood location of longer-duration energy-release source regions, active during the tectonic tremor sequence (Fig. 10a).
The broad-band envelope CFs might however still contain remaining short-scale fluctuations, possibly originating from impulsive short-duration low-frequency earthquakes within the 180 s sliding time-windows (see Supporting Information Fig. S4a).The LCC-based TDE computation (Section 4.1.2),introduces a 2-D Gaussian filter with a half-width σ , in the time-lag space, which defines the coherency scales dominating the LCC.For the case of longer-duration energy-release source analysis this Gaussian filter needs to be adapted to emphasize and extract within the sliding time-widows longer-duration coherency scales between stationpairs CFs, smoothing out the influence of the short-scale fluctuations (Supporting Information Figs S4b, d and e).This procedure might be related to the low-pass filtering of the energy envelopes used in other studies (e.g.Ide 2010).The influence of LCC 2-D Gaussian filtering on station-pair TDE distributions, the imaging function and the retrieved final source location are illustrated in Supporting Information Fig. S4.The example assumes different values of σ , that is, σ = 0 s and σ = 3.5 s.In the case of σ = 3.5 s, the LCC function is dominated by long-duration fluctuations of the CFs, which provides a 3-D likelihood imaging function corresponding to an extended energy-release source region.In this study, a fixed value of σ = 3.5 s is selected allowing to focus on the detection and imaging of longer-duration energy-release sources during the analysed 9-day tectonic tremor sequence.
Once the imaging function for a given time-window is estimated, the detection of a potential longer-duration energy-release source is based on two criteria: the first one involves the value of the maximum likelihood of the imaging function, while the second one combines the maximum likelihood and the size of the major axis of the 68 per cent confidence ellipsoid (trigger ellipsoid), estimated with respect to the (x, y, z) mean point of the imaging function (Lomax 2005(Lomax , 2008)).Accordingly, detection of a potential longer-duration energy-release source is declared when the maximum likelihood value of the imaging function is ≥0.95 (95 per cent) or when the major axis of trigger ellipsoid is ≤50 km and the maximum likelihood value of the imaging function is ≥0.85 (85 per cent; Fig. 10a).This allows detecting low-energy sourcesrecorded by few stations-when those correspond to a well-focused imaging function, while decreasing the number of false detections.In contrast to low-frequency earthquakes analysis-for which the point-source location of an impulsive short-duration event is associated with the maximum likelihood point in the imaging functionhere, the entire imaging function is retained.We, therefore, refer to the procedure as detection and imaging, rather than detection and location.The output of this step are the time-series of imaging functions that provide space-and-time likelihood information on the location of longer-duration energy-release sources during the tectonic tremor sequence, that is, space-and-time evolution of the energy-release region for a given likelihood threshold.The imaging functions provide a smoothed representation of the effective longer-duration source region that includes also the arrayresponse function of the seismic network.
The final phase of the detection and imaging scheme involves a post-processing step that is applied to the time-series of imaging functions in order to remove redundant information, associated with the overlap between consecutive sliding time-windows, and to account for the spatial and temporal extension of the longerduration energy-release sources.This bares some similarity with the underlying idea of the clustering scheme proposed by Obara et al. (2010), and implemented in the NIED tectonic tremor catalogue.The post-processing consists in associating and stacking the imaging functions that correspond to detections in at least 3 successive sliding time-windows, and to horizontal locations of the maximum likelihood that are no more than 15 km apart while the origin times are within 20 min difference, see Fig. 10 for illustration.An example of such summary imaging function corresponding to a single longer-duration event is given Fig. 10(b).This final step leads to a catalogue of longer-duration energy-release events characterized by their time-duration and the likelihood information on the energy-release source region described by the corresponding summary imaging function.These summary imaging functions can be further combined in time, that is, stacked, to form a 9-day summary imaging function that provides the likelihood of longer-duration energy-release sources in a given region during the period of tec-tonic tremor sequence (Fig. 11), or which can be analysed in terms of the temporal evolution of the longer-duration energy-release sources.Alternatively, a reduced and more classical catalogue of effective point-source locations can also be obtained by extracting the spatial position and the time where the maximum likelihood is achieved in the summary imaging functions.The time-series of summary imaging functions provide rich space and time information about longer-duration energy-release source regions, during tectonic tremor sequence, in the frequency range of 2.0-15.0Hz, as a result of the broad-band energy CFs reconstruction (Section 5.1.1,Fig. 9).

Results
The 9-day summary imaging function, covering the entire period of the tectonic tremor sequence, is shown in Fig. 11.The largescale longer-duration energy-release source region-during the 9-day tectonic tremor sequence-is characterized by a predominant southwest-northeast extension, of about 150 km, along the strike direction of the subducting slab.The extension in the direction perpendicular to the strike is of about 50-75 km.Even though depth-dependent information is less constrained, most of the emitting source region appears to be localized in zones 2 and 3 in Fig. 11, concentrated in depth between 10-60 km.This observation is relatively less evident in the southwestern part of the western Shikoku, which might be related to a more spread activity at the beginning of the tectonic tremor sequence previously evidenced by the low-frequency earthquakes catalogue.Nonetheless, one of the main characteristics of the 9-day summary imaging function is a clear along-strike segmentation of the longer-duration energy-release source regions.While the summary imaging-function is relatively diffuse-with no clear likelihood maxima-in the southwestern part (zone 1 in Fig. 11), the central and north-eastern parts-zones 2 and 3 in-exhibit well defined patches of high energy-release activity.
The setup of the BackTrackBB method applied here is designed to extract and restore longer-scale spatial and temporal features of energy-release during the 9-day tectonic tremor sequence.As such, slow along-strike migration of the energy-release activity could be recovered from the time-series of summary imaging functions, while shorter and smaller-scale migration patterns can only be observed from the low-frequency earthquakes catalogue.
Comparison between the 9-day summary imaging function and the tectonic tremor catalogue of NIED, providing 1 hr averaged effective point-source locations of long-duration energy sources during the tectonic tremor sequence is shown in Fig. 11.The location and the extension of the high energy-release patches extracted from the 9-day summary imaging-function (e.g.zones 2 and 3 in Fig. 11) are in good agreement with the clusters of tectonic tremor point source epicentres from the NIED catalogue.Most of the NIED epicentres are located well within the energy-release regions of maximum likelihood of the summary imaging function.The distribution of the NIED epicentres is relatively spread towards the southwest in agreement with the more diffuse appearance of the summary imaging function in the same region (zone 1 in Fig. 11).This suggests that the summary imaging functions do capture and resolve large-scale properties of the longer-duration energy-release sources during the tectonic tremor sequence despite the poor station coverage in this region.
The 9-day summary imaging-function can also be analysed together with the catalogue of low-frequency earthquakes derived in this study (Fig. 12a).The spatial extension and the segmentation of the long-duration energy-release source regions are consistent with observations derived from the analysis of the low-frequency earthquakes (Figs 4a and 5a).In particular, the high energyrelease patches extracted from the 9-day summary imaging function (zones 2 and 3 in Fig. 11) agree well with the clusters of lowfrequency earthquake activity (Figs 4a and b) leading to the conclusion that both short-scale impulsive low-frequency earthquakes and longer-scale energy-release sources, during the 9-day tectonic tremor sequence, appear to occupy the same regions.Moreover, high-activity regions of low-frequency earthquakes are associated with the high energy-release zones of longer-duration energy sources, and large-scale spatial and temporal patterns of both longer-duration source activity and low-frequency earthquake activity regions are well correlated.To understand better the relationship between the short-duration (low-frequency earthquakes) and longer-duration source components during the tectonic tremor sequence, we carried out a simple preliminary statistical analysis of the corresponding catalogues.We calculated and examined the distribution of distances between the centroids of the summary imaging functions associated with individual longer-duration energy-release events (Fig. 10b) and the centroids of low-frequency earthquakes detected during these events.We looked as well at the distances between the individual low-frequency earthquakes and the longer-duration events' centroids.The results of these investigations are summarized in the histograms presented in Figs 12(b) and (c).More details regarding the statistics of the number of detected lowfrequency earthquakes and the duration of extracted longer-duration events can be found the Supporting Information (Fig. S5).There are in total 313 longer-duration events (and their corresponding imaging functions) identified during the tectonic tremor sequence in our final catalogue of longer-duration sources.Forty-eight of these events do not have corresponding low-frequency earthquakes detections, which raises the question of the actual detectability of low-frequency earthquakes during these periods or of the origin of this energy-release component.From 2212 low-frequency earthquakes in the final catalogue, 1655 are appearing during the longerduration events.The investigation of the histogram in Fig. 12(b) indicates that the low-frequency earthquakes centroids and the corresponding longer-duration event centroids are coincident within about 10-15 km distance.This estimate corresponds to an upper limit, due to the respective location uncertainty of the lowfrequency earthquake and longer-duration energy-release sources.The result is stable when looking at the distances to the individual low-frequency earthquakes (Fig. 12c) and is also supported by histogram in Fig. 12(b), indicating that most of the low-frequency earthquakes fall inside the region of summary imaging-function >0.98 (Fig. 12d).These suggest a physical correlation, potentially involving short-scale interactions, between the short-duration (low-frequency earthquakes) and the longer-duration energy-release source components during the tectonic tremor sequence.Moreover, low-frequency earthquakes are observed as localized in the vicinity of the longer-duration energy-release events centroids implying that the highest energy-release extended sources during tremor activity are likely to be restricted in size to a radius of less than 10-15 km.This provides a first quantitative estimate for the size of the longer-duration energy-release sources, such as, tremor sources of the NIED catalogue, during a tectonic tremor sequence in western Shikoku.

D I S C U S S I O N S A N D C O N C L U S I O N S
In this study, analysing a 9-day long energetic tectonic tremor sequence, we showed that extracted short-duration low-frequency earthquakes activity is coincident in space (within about 10-15 km distance) with the longer-duration energy-release source component.Spatial location of tectonic tremor is an open problem (Rubinstein et al. 2009;Ide 2012) due to the possible extended size and heterogeneity of the emitting source area that might be scale bound (Watanabe et al. 2007) rather than scale invariant.While lowfrequency earthquakes may contribute to self-excitation of tremors, longer-duration energy components in the signal of a tectonic tremor sequence can possibly represent the signature of nonlinear mixing of energy-release sources that remains largely unknown today.Our analysis pointed out that both impulsive short-duration and longerduration source components during this tectonic tremor sequence occupy the same (i.e.coincident within 10-15 km distance) effective energy-release source regions (Figs 12a-c), providing a quantitative constrain on the size of the energy radiating source areas during the tectonic tremor sequence.This result, however, inferred from the case of one tectonic tremor sequence in the Nankai subduction zone, might not be relevant for the ETS activity in other subduction regions, such as Mexico or Cascadia.
In this study, we applied an extended formulation of the recently developed network-based BackTrackBB method (Poiata et al. 2016) to the analysis of a 9-day long (from 2012 May 25 to June 2) tectonic tremor sequence.Through the analysis presented here, we illustrated new capabilities of the method to automatically detect and locate different scale components of the energy-release activity during an energetic tectonic tremor sequence in southwestern Japan.We performed an automatic detection and location of the shortduration, impulsive low-frequency earthquakes using characteristic functions (CFs) based on frequency-dependent HOS of the signal recorded at each station of the network.The obtained automatic BackTrackBB catalogue of low-frequency earthquakes provides a more comprehensive and time-resolved information on the lowfrequency earthquake activity during the analysed tectonic tremor sequence than the existing JMA catalogue of manually detected low-frequency earthquakes.
The automatic BackTrackBB catalogue allows imaging the space-time activity of the low-frequency earthquakes during the tectonic tremor sequence in great detail, unravelling spatial and temporal clustering, modulation of activity in response to tide, and different scales of space-time migration patterns: a large-scale along-strike slow propagation of the main front of low-frequency earthquakes with an average speed of 10-12 km d −1 ; a shorterscale and faster reverse-strike migration behind the main front with speeds of 5-7 km hr −1 ; and a short-scale and rapid along-dip migration with speeds varying from 7-10 to 75-80 km hr −1 .Slow along-strike propagation of the low-frequency earthquakes activity front is well known in the Shikoku area and other regions worldwide (e.g.Obara 2002;Dragert et al. 2004), and previous studies (Shelly et al. 2007a,b) already evidenced fast along-dip migration patterns in Shikoku from template match-filtering techniques.To our knowledge, this is the first time that reverse-strike migration of the low-frequency earthquakes activity, previously observed by Huston et al. (2011) in Cascadia, is evidenced in this area from a single 9-day tectonic tremor sequence.It is also worth noting here that, in contrast to match-filtering technique, BackTrackBB does not require pre-existing template dictionary and does not assume a priori that low-frequency earthquakes activity during a tectonic tremor sequence is fully characterized by families of repeated sources.As such, BackTrackBB is an attractive automatic high-resolution detection and location method that might complement the match-filtering approach providing efficient mean to construct comprehensive template dictionaries.Some details about the computational aspects and efficiency of BackTrackBB in regard to the presented analysis are presented in Supporting Information Section S1.
In the second part of the study, longer-duration components of low signal-to-noise ratio energy transients during the tectonic tremor sequence were analysed using CFs built from time-frequency energy envelopes.BackTrackBB allows to efficiently detect longerduration energy transients and to image their extended energyrelease source regions.The resulting catalogue consists of timeseries of summary imaging functions providing 3-D likelihood maps of longer-duration energy-release events together with their durations.This is in contrast to classical methods, and existing catalogues, providing only averaged effective tectonic tremor pointsource locations (e.g.Maeda & Obara 2009;Ide 2010).
The 9-day summary imaging function -covering the entire period of the tectonic tremor sequence-provides information about 3-D spatial distribution of extended long-duration energyrelease regions during the sequence, which are characterized by clear along-strike segmentation.Comparison with the low-frequency earthquakes activity shows that, along the strike of subducting plate, high energy-release zones of longer-duration energy sources largely coincide with clusters of intense low-frequency earthquakes activity.This clustering and segmentation might, as proposed in studies (e.g.Shelly et al. 2007a;Ide 2010), be associated with 'frozen' heterogeneity along the plate boundary as observed in other subduction zones such as Cascadia (e.g.Ghosh et al. 2010b) and Guerrero, Mexico (e.g.Peng & Rubin 2017) The large-scale along strike migration patterns of both longduration energy-release source areas and low-frequency earthquakes activity are actually quite similar during the 9-days of tremor sequence.However, smaller-scale spatial and temporal features can be only resolved from the analysis of the low-frequency earthquakes catalogue that provides high-resolution information of short-scale energy-release processes.This is well exemplified by the downdip extension of low-frequency earthquakes hypocentres around the centre of the active long-duration energy-release source region (Figs 6 and 12a), and the low-frequency earthquake cluster located between zones 1 and 2 evidenced in the summary imaging function (Fig. 11).It is yet unclear whether this might be related to the difference in scale resolution between low-frequency earthquakes and the long-duration energy-release source analysis, or it is evidencing differences in short-scale and longer-scale dynamics during the tectonic tremor sequence.Some of the emitting long-duration energy-release source regions might also be more prone to triggering low-frequency earthquakes.Future systematic analysis applying the BackTrackBB method to larger data sets, containing multiple tectonic tremor sequences, might answer this question.
The construction of the energy envelope CFs using the MBF decomposition (Section 5.1.1 and Fig. 9) opens also new interesting directions for investigating frequency-dependent components of tectonic tremor signals.This can be achieved introducing into the detection and location step CFs that are built from only selected sub-bands of frequency-dependent energy envelopes after the MBF signal decomposition.An illustration of such a potential application is provided in Supporting Information Fig. S6.The figure shows the imaging functions and the summary imaging function associated with a detected tremor energy-release source for the selected MBF sub-bands: 2.00-2.68Hz (Supporting Information Fig. S6b), 6.11-6.78Hz (Supporting Information Fig. S6c) and 10.21-10.89Hz (Supporting Information Fig. S6d).It can be noted that the predominant frequency component, defining the scales of the energy envelope fluctuations, is reproduced, as well, in the imaging functions.Such frequency-dependent analysis of the tremor energy-release sources provides a way to image different scales of the tectonic tremor energy-release regions both in space and time, and might lead to original insights into potential scale-dependent dynamics of the underlying processes.This is, however, out of the scope of this paper, remaining as the suggested possible extension of this work.
One of the conclusions of the present study, that we would like to emphasize, is that network-based methods using advanced statistical signal processing and coherent imaging schemes, such as BackTrackBB, offer today new capabilities for high-resolution detection, location and analysis of the different scale-components of tectonic tremor activity, enriching existing slow earthquakes catalogues.The data streaming BackTrackBB method offer also new capabilities for continuous monitoring of tectonic tremor activity.Application of the method to the continuous data set, covering years of recording in southwestern Japan, might lead to new observations and insights into the underlying scale-dependent mechanisms of slow earthquakes.

Figure 1 .
Figure 1.(a) Map view of the southwestern Japan study-area; grey and black triangles indicate location of Hi-net (NIED) seismic stations.Grey circle corresponds to the location of the tectonic tremor source reported by NIED for the 1 hr period of 0300-0400 JST, 2012 May 31.Yellow star corresponds to the epicentre of a selected low-frequency earthquake (t 0 = 03:19:31.89,2012 May 31) according to JMA catalogue.(b) Record section of 1 hr NS component seismograms for the tremor activity corresponding to the tremor source epicentre shown in (a).The seismograms from selected stations-black triangles in (a)-are filtered with a bandpass of 2.0-25.0Hz, normalized and ordered according to the increasing epicentral distance (scale is not preserved, and added only as reference) from tremor event source in (a).Red line corresponds to a selected low-frequency earthquake.(c) Expanded view of the seismogram showing a low-frequency earthquake during tectonic tremor activity in (b) reported by JMA catalogue.Red lines indicate the theoretical S-waves arrivals calculated according to the JMA epicentre location (a), and the velocity model of Kubo et al. (2002).

Figure 2 .
Figure 2. Illustration of MBF statistical signal processing scheme applied for estimating the characteristic functions that extract and enhance information on low-frequency earthquakes.(a) Filtered traces.(b) HOS kurtosis CF of filtered signals.(c) Raw broad-band seismogram.(d) Broad-band HOS kurtosis CF (in red) composed as the maxima over all frequencies, and the final broad-band MBF CF (in blue), obtained by convolving the CF's positive derivative with a Gaussian window.(e) Comparison of HOS Kurtosis CF (black) and the final CF (dotted black) calculated on raw data, and the broad-band MBF CF's (red and blue).Grey shadowed window in (c-e) corresponds to the selected low-frequency earthquake shown in Fig. 1(c).

Figure 3 .
Figure 3. Illustration of the detection and location scheme for one example of selected low-frequency earthquake.(a) Final broad-band MBF CFs (blue lines) of a selected time-window containing low-frequency earthquake (same as Fig. 1c and Fig. 2).The NS components of the original velocity seismograms are shown by grey lines.Traces are arranged from bottom to top by the order of increasing epicentral distance (distance scale is not respected).Vertical red bars indicate theoretical S-wave arrival times from estimated location.(b) One horizontal and two vertical sections through the final estimated location, corresponding to the maximum of imaging function, represented using a colour scale.White triangles indicate the location of Hi-net stations used for location; green star indicate the resulting location of low-frequency earthquake.(c) Zoomed image of panel (b) comparing the imaging function and the location of the event provided by JMA catalogue (white star).

Figure 4 .
Figure 4. (a) Map of western Shikoku showing the distribution of low-frequency earthquakes (coloured circles) from the catalogue obtained by automatic analysis of 9-day data set.Circle colour corresponds to the timing of the events relative to the beginning of the tectonic tremor sequence on 2012 May 25, as indicated by the colour scale.Black squares show the location of the Hi-net stations used for analysis.Dashed lines are the depth contours of the Philippine Sea Plate top (PSt) surface derived from the Japan Integrated Velocity Structure Model (Koketsu et al. 2008).(b) Along-strike (A-A , N35 • E) cross-section of the detected low-frequency earthquakes as the function of depth.Circle colour corresponds to colour scale in (a).Dashed line represents the depth of the Philippine Sea Plate top surface along the cross-section A-A , from Koketsu et al. (2008).Solid black line shows the average depth profile fitted to the along-strike 4 km binned low-frequency earthquakes hypocentres.(c) Along-dip cross-sections of the detected low-frequency earthquakes corresponding to the different selected profiles in (a).Dashed lines correspond (from left to right) to the mantle wedge limit (MW), top and bottom of the Philippine Sea Plate (PSt and PSb respectively), after Koketsu et al. (2008).

Figure 5 .
Figure 5. Side by side comparison between the automatic low-frequency earthquake catalogue derived in this study and the JMA low-frequency earthquake catalogue.(a) Distribution of the low-frequency earthquakes (coloured circles) from the catalogue obtained by automatic analysis of the 9-day tectonic tremor sequence shown on the map of western Shikoku and as two depth-dependent cross-sections: along-strike (A-A ) and along-dip (perpendicular to A-A ).(b) Hourly activity rate of the detected low-frequency earthquakes during the analysed 9-day tectonic tremor sequence, in number of events per hour.Histogram is colour-coded according to the colour-scale shown below the figure.(c) Distribution of the low-frequency earthquakes reported by the JMA catalogue for the same 9-day time-period, shown on the map of western Shikoku and two depth-dependent cross-sections: along-strike (A-A ) and along-dip (perpendicular to A-A ).(d) Hourly activity rate of low-frequency earthquakes reported by the JMA catalogue, in number of events per hour.

Figure 6 .
Figure 6.(a) Space-time plot of the detected low-frequency earthquake activity, showing the along-strike (A-A ; Fig. 5a) distribution of events with time.Circle colour corresponds to the timing of the low-frequency earthquakes, relative to the beginning of the tectonic tremor sequence on 2012 May 25, as indicated by the colour scale.Histogram in grey corresponds to hourly activity rate of detected low-frequency earthquakes in number of events per hour.Black arrows indicate the directions of along strike main front migration and reverse-strike migrations during the intensive bursts of activity marked by grey rectangles that are numbered consecutively.(b) Space-time plot of the low-frequency earthquake activity showing the along-dip (perpendicular to A-A ; Fig. 5a) distribution of events with time.Histogram in grey corresponds to hourly activity rate of detected low-frequency earthquakes in number of events per hour.(c) Normalized 9-day energy envelopes of the NS component velocity seismograms recorded at three Hi-net stations.

Figure 7 .
Figure 7. Map view and space-time plots of the low-frequency earthquakes extracted from the automatic analysis and covering a 23 hr selected time-window of intense activity (period number 3 in Figs 5 and 6).(a) Map of the low-frequency earthquake epicentres for the selected 23 hr period (coloured circles).Circle colour corresponds to the timing of the low-frequency earthquakes relative to the beginning of the selected time-period on 2012 May 30, as indicated by the colour scale.Light grey circles correspond to the entire catalogue of low-frequency earthquakes detected and located in this study.(b) Hourly activity rate of the detected low-frequency earthquakes during the selected period, presented as number of events per hour.Colour scale is same as in (a).(c) Cumulative number of detected low-frequency earthquakes.(d) Along-strike (A-A ; Fig. 5a) distribution of the low-frequency earthquakes with time during the selected period.(e) Along-dip distribution of the low-frequency earthquakes as a function of time during the selected period.

Figure 8 .
Figure 8. Normalized 23 hr NS component velocity seismograms recorded at three Hi-net stations.The red highlights indicate the time-windows corresponding to the low-frequency earthquakes located in this study and shown in Fig. 7.

Figure 9 .
Figure 9. Illustration of the MBF signal processing applied for the estimation of the energy envelope CFs.(a) Filtered traces.(b) Energy envelope CFs of filtered signals.(c) Raw broad-band seismogram.(d) Broad-band CF (in red) composed with RMS operator applied over all the frequencies of MBF filter, and 2-15 Hz filtered seismogram (in grey).

Figure 10 .
Figure 10.(a) Imaging functions for three selected time-windows, part of the consecutive time-window for which the long-duration energy source was detected and associated to a single event.The black triangles show the Hi-net stations used in the analysis.White ellipses correspond to 68 per cent confidence ellipsoid (trigger ellipsoid), estimated with respect to the (x, y, z) mean point of the imaging function.(b) Summary imaging function estimated as the normalized sum of the imaging functions for multiple time-windows (like those in panel a), associated as corresponding to a single long-duration energy-release event.White circle indicates the maximum of summary imaging function.Small yellow stars correspond to the location of low-frequency earthquakes detected by the automatic scheme within the same time-window.(c) Normalized 1 hr NS component velocity seismograms recorded at three Hi-net stations (black lines) and corresponding broad-band envelope characteristic functions (red lines).Green shaded rectangle corresponds to the time-windows with consecutive detections of long-duration sources associated to the summary imaging function in (b).Rectangles delimited by thick green lines indicate the three selected time-windows corresponding to the imaging functions shown in (a).

Figure 11 .
Figure 11.9-day summary imaging function of the detected longer-duration energy-release sources during the tectonic tremor sequence.Black triangles indicate Hi-net stations.Grey open circles indicate the location of the tectonic tremor point sources reported for the same sequence in NIED tremor catalogue.

Figure 12 .
Figure 12.(a) Map and two vertical cross-sections comparing the 9-day summary imaging function (SLF) and the distribution of the low-frequency earthquakes from the catalogue obtained by automatic analysis of 9-day data set in this study (grey dots).Black triangles show the location of Hi-net stations.(b) Histograms of the horizontal distances between the centroid location of the detected longer-duration events (LDEs) and the centroid of the low-frequency earthquakes (LFEs) occurring during the longer-duration events within the 9-day tectonic tremor sequence.(c) Histograms of the horizontal distances between the centroid location of the detected LDEs and the individual LFEs occurring during these events within the 9-day tectonic tremor episode.(d) Histogram indicating the LDEs value of summary imaging function (SIF) at the location of low-frequency earthquakes occurred during the duration of the event.