Cascading rupture of patches of high seismic energy release controls the growth process of episodic tremor and slip events

Slip phenomena on plate interfaces reflect the heterogeneous physical properties of the slip plane and, thus, exhibit a wide variety of slip velocities and rupture propagation behaviors. Recent findings on slow earthquakes reveal similarities and differences between slow and regular earthquakes. Episodic tremor and slip (ETS) events, a type of slow earthquake widely observed in subduction zones, likewise show diverse activity. We investigated the growth of 17 ETS events beneath the Kii Peninsula in the Nankai subduction zone, Japan. Analyses of waveform data recorded by a seismic array enabled us to locate tremor hypocenters and estimate the migration patterns and spatial distribution of the energy release of tremor events. Here, we describe three major features in the growth of ETS events. First, independent of their start point and migration pattern, ETS events exhibit patches of high seismic energy release on the up-dip part of the ETS zone, suggesting that the location of these patches is controlled by inherent physical or frictional properties of the plate interface. Second, ETS events usually start outside the high-energy patches, and their final extent depends on whether the patches participate in the rupture. Third, we recognize no size dependence in the initiation phase of ETS events of different sizes with comparable start points. These features demonstrate that the cascading rupture of high-energy patches governs the growth of ETS events, just as the cascading rupture of asperities governs the growth of regular earthquakes.


Introduction
Earthquakes feature heterogeneous slip distributions on fault planes that are characterized by one to several areas of large slip, or asperities. In general, earthquake hypocenters are located outside of asperities (e.g., Yamanaka and Kikuchi 2004;Mai et al. 2005). The complex growth and rupture processes of earthquakes are fundamental features that have fueled a long-running debate on how and when the eventual size of an earthquake is determined, specifically the existence of a nucleation phase (e.g., Ellsworth and Beroza 1995;Ide 2019).
Recent findings about slow earthquakes have broadened our understanding of subduction dynamics and the generation of megathrust earthquakes (e.g., Obara and Kato 2016). In subduction zones, strain on the plate interface is released by slow slips in slow earthquakes and high-speed ruptures in regular earthquakes. Interestingly, the spatial distributions of slow and regular earthquakes are complementary: slow earthquakes are distributed around and outside the locked zones that are the source areas of megathrust earthquakes, both along dip (Obara and Kato 2016) and along strike (Nishikawa et al. 2019). Furthermore, a temporal relationship between slow and regular earthquakes has been proposed. Kato et al. (2012) reported the occurrence of two slow slip events (SSEs) from earthquake activities that migrated toward the future hypocenter of the 2011 Tohoku earthquake before the mainshock. The migration speed of those hypocenters, 2-5 km/day, is close to the typical migration speed of episodic tremor and slip (ETS) events, 10 km/day (Dragert et al. 2001;Obara 2002). Slow slips preceding a dynamic rupture have also been reported in laboratory experiments and numerical simulations (e.g., Dieterich 1979;Ohnaka 1992;Shibazaki and Matsu'ura 1992). Matsuzawa et al. (2010) showed that in a numerical simulation, the recurrence interval of SSEs grew shorter in the time interval preceding a megathrust earthquake. Therefore, investigating the features of slow earthquakes yields clues to the generation of megathrust earthquakes in subduction zones.
In some subduction zones, slow earthquakes occur as short-term SSEs in the form of ETS events, in which tectonic (nonvolcanic) tremor and low-frequency earthquake (LFE) are synchronized with slip events in time and space (Rogers and Dragert 2003). The slip areas and slip rates of ETS events are strongly correlated with tremor epicenters Bartlow et al. 2011;Michel et al. 2019a). Thus, tremor and LFE can be used to monitor ETS events (Shelly et al. 2007a;Obara 2010;Houston et al. 2011;Ghosh et al. 2012;Rubin and Armbruster 2013;Peng et al. 2015), although short-term SSEs are not always accompanied by tremors (Obara et al. 2011;Wech and Bartlow 2014;Michel et al. 2019a). Tremor and LFE amplitudes are also have been used to detect small episodes of transient aseismic slip hidden in the geodetic noise (Frank 2016;Frank et al. 2018). The growth process of ETS events is not simple and is typified by diverse patterns of migration and eventual sizes , heterogeneous slip distributions , and brief episodes of secondary slip front in which tremor migrates faster than the main front: sometimes in the direction opposite to its overall migration (Houston et al. 2011). Ghosh et al. (2012) reported the existence of tremor asperities, areas of high radiated seismic energy, on the plate interface analogous to the asperities of regular earthquakes and suggested that tremor asperities in the transition zone might control the evolution of slow earthquakes in space and time. Savard and Bostock (2015) identified the location of five asperities from LFE epicenters in Cascadia. Hirose and Obara (2010) conducted geodetic inversion and confirmed the existence of a persistent asperity of slow slip approximately 30 km in diameter in Shikoku, Japan, that was active during several different ETS events. In the Nankai subduction zone off southwestern Japan, the ETS zone can be divided into several segments (e.g., Obara 2002) within which ETS events have similar recurrence intervals, such as 3 or 6 months. The short recurrence intervals of these ETS events make the Nankai subduction zone a promising area to elucidate what controls their growth.
The moment rate of slow slips is often estimated from the number, amplitude, or radiated energy of tremors or LFEs. Ide et al. (2008) demonstrated that the radiated seismic energy rate of tremors in the 2-8 Hz band is proportional to the seismic moment rate of very-lowfrequency earthquakes in the Kii Peninsula. Ide and Yabe (2014) confirmed this proportionality between tremor energy rate and seismic moment rate of very-low-frequency earthquakes in the entire Nankai subduction zone. Hirose and Obara (2010) showed that the moment rate function of an SSE estimated from tiltmeter records was strongly correlated with the temporal change in the number of tremors in Shikoku. Maeda and Obara (2009) found a clear correspondence between the temporal variation in tilt records and the accumulated seismic energy from tremor sources for ETS events in Shikoku. Hawthorne and Rubin (2013) found that the strain rate is proportional to the amplitude of tremor on time scales shorter than one day. They interpreted this relation to mean that the moment rate of slow slip is correlated with tremor amplitude.
Seismic arrays are a powerful tool for detecting, locating, and characterizing tremor (Ghosh et al. 2009;Imanishi et al. 2011). The Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology (GSJ), installed a seismic array of 39 sensors on the Kii Peninsula from February 2011 to November 2016 ( Fig. 1), producing a set of well-recorded waveforms that has been valuable for examining the detailed growth process of ETS events from tectonic tremor. Sagae et al. (2021) applied the MUSIC high-resolution frequency-wavenumber (f-k) method (Schmidt 1986) to a dataset from July 2012 to July 2014, revealing a more detailed view of tremor migration than the conventional envelope crosscorrelation method. We report here the spatio-temporal distribution of the radiated seismic energy of tectonic tremor during 17 ETS events between April 2011 and December 2014. We show that persistent patches of relatively high seismic energy release are located on the plate interface in the ETS zone, and that the growth process of ETS events is controlled by cascading ruptures of these patches.

Data
We analyzed waveform data recorded by the GSJ array network in the Kii Peninsula from April 2011 to December 2014 (Fig. 1). The array consisted of 39 three-component velocity seismographs with a natural frequency of 2 Hz, and the sampling frequency of waveform data was 200 Hz. For semblance analysis, we used the N-S component of the seismic data bandpass-filtered between 2 and 4 Hz because the horizontal components yielded higher semblance value than the vertical component. For the calculation of the seismic energy of tremors, we used all three components of data between 2 and 8 Hz. We used the JMA2001 seismic velocity model (Ueno et al. 2002) to locate tremors.

ETS event durations for semblance analysis
In general, the size and temporal evolution of an ETS event are correlated with the number of recorded tremors (e.g., Obara 2010). In this study, we treated a concentration of tremors in space and time as an ETS event. Tremor locations were determined by semblance analysis (Neidell and Taner 1971), in which we estimated the horizontal slowness vector of coherent signals under the assumption of a plane-wave incidence, as described in the next subsection. To estimate the duration of ETS events from tremor activity for our semblance analysis, we applied the following clustering technique to the National Research Institute for Earth Science and Disaster Resilience (NIED) hybrid clustering catalog , which locates one tremor per hour. We classified tremor events as hypocenters clustered within 20 km during a 48-h period. Groups of more than 10 hypocenters were considered to be ETS events. We estimated  (black rectangle;Figs. 5,6,7,8) and the array network (red square). Gray dots are tremor hypocenters from the NIED hybrid clustering catalog. Contour lines represent the depth of the Moho in the subducting Philippine Sea plate (Shiomi et al. 2008). Dotted square represents the area shown in (c). b Configuration of stations (triangles) in the seismic array. The black triangle is the central station. c The distribution of short-term SSEs detected geodetically by GSJ and/or NIED during the study period. Gray rectangles are surface projections of areas on the fault planes of the short-term SSEs; the black side of each rectangle is the up-dip edge. Yellow triangle represents a Hi-net station (KRTH) for which tilt data are shown in Fig. 4b. d Space-time plot of tectonic tremors from the NIED hybrid clustering catalog. Distance along strike or along dip is measured from the corner with a black circle of the black rectangle in (a). Circles represent hypocenters of tremor clusters for different ETS events analyzed in this study, colored according to their ETS number. ETS events 1, 8, and 11 were excluded from further analysis. Gray dots are tremor hypocenters that are isolated or belong to small ETS events. Gray vertical bars represent short-term SSEs detected geodetically by GSJ the duration of an ETS event as the time from the start of the group's first tremor to the end of its last tremor. Consequently, we detected 20 ETS events from June 2011 to July 2014 for our semblance analysis, of which 16 had previously been geodetically detected as short-term SSEs from analyses of strainmeter, tiltmeter, and/or water-level records Hirose and Kimura 2012;Itaba et al. , 2013aItaba et al. , 2013bItaba et al. , 2014aItaba et al. , 2014bItaba et al. , 2015Ochi et al. 2015) (Figs. 1c and 1d). The semblance analysis utilized array data from two days before the start to one day after the end of each ETS event (Table 1).

Locating tremors from semblance analysis
In this study, we used waveform data recorded by the GSJ array to detect more tremors. Tremor has been shown to occur mainly on the plate interface (e.g., Shelly et al. 2007b;Ghosh et al. 2012;Suzuki et al. 2018). Thus, we assigned possible hypocenters to grid points with horizontal intervals of 2 km on the plate interface, which is 5 km shallower than the Moho in the subducting Philippine Sea plate (Shiomi et al. 2008). We calculated the slowness of the optimum ray path between each grid point and the central station of the array network within the JMA2001 velocity model (Ueno et al. 2002).
Furthermore, we compensated for the arrival time shift due to the elevation difference between the stations. We used eight earthquakes with magnitudes ≥ 1.0 and epicentral distances from the central station of the array network ≤ 10 km from 2011 to 2017 (Fig. 2a). We calculated the arrival time differences between the central station and the other stations from the cross-correlation function of waveform data using a time window of 1 s around the P-wave arrival time. The arrival time differences are correlated positively with the elevation differences for all earthquakes. We assumed here that P-wave ray paths were vertical directly beneath the array network, so that the arrival time differences depended only on the elevation differences. A linear regression line between the  arrival time differences and the elevation differences was determined by the least-squares method constrained to pass through the origin (Fig. 2b). The slopes of the lines, as P-wave velocities, were highly variable because of a large scatter in the data (Additional file 1: Fig. S1). Therefore, we chose the data of the earthquake that showed the best fitting (Fig. 2b). As a result, the estimated P-wave velocity was 3.0 km/s (Fig. 2b), providing a S-wave velocity of 1.7 km/s with the assumption of V P /V S = √ 3 . We applied this S-wave velocity above sea level for the semblance analysis. We confirmed that the variation (± 30%) in this S-wave velocity did not significantly modify the following results and discussion.
We conducted semblance analysis (Neidell and Taner 1971) of the continuous records to detect and locate tremors. Compared to the conventional beamforming and f-k methods, semblance analyses can detect events with lower signal-to-noise ratios because the method is based on waveform coherency rather than the signal power (see the review by Rost and Thomas 2002, for the advantages and disadvantages of various array techniques). From the waveform data of the N-S component, bandpass-filtered at 2-4 Hz, we calculated semblance from the slowness values of the grid cells by where SEM is a semblance value, W a time window of 1 min, T the start time of a time window, M the number of stations for which the maximum cross-correlation value of waveforms for time shifts of -0.3 to + 0.3 s between the central station and the j-th station is greater than 0.3, f j the waveform data at the j-th station, τ nj the time difference between the central station and the j-th station calculated from the horizontal slowness vector (S nx , S ny ) of the n-th grid point, and x j and y j are the x and y coordinates of the j-th station, respectively. HC nj is the time correction for the elevation difference between the central station and the j-th station for the horizontal slowness vector of the n-th grid point, and SC j is the station correction that is estimated to match hypocenters from the catalog with those from the semblance analysis.
If the maximum semblance value was greater than 0.3 and the number of stations used for the analysis was greater than 25, the grid point with the maximum semblance value was designated the hypocenter of a tremor. We estimated the uncertainty of tremor locations by applying the delete-1 jackknife method (Ueno et al. 2010). The estimated uncertainty depends on the epicentral (1) (2) τ nj = S nx · x j + S ny · y j − HC nj − SC j , distance from the array; the typical uncertainty is 2-4 km for epicentral distances less than 35 km, 4-6 km for distances of 35-50 km, and 6-10 km and more for distances greater than 50 km. We also recognize that tremor locations estimated by the delete-1 jackknife method tend to be distributed in a region that is elongate along the backazimuth direction at epicentral distances greater than 40-50 km.
By applying a clustering procedure to the estimated hypocenters, we removed isolated events that were not part of an ETS event. Our criteria for clustering two hypocenters were an epicentral distance less than 5 km and a time interval shorter than 12 h. We visually inspected the waveforms of detected tremors to exclude waveforms of regular earthquakes or impulsive noises.

Estimation of the seismic energy of tremors
We calculated seismic energy from waveforms of three components at the central station during 1-min time windows by following the formula (Maeda and Obara 2009) where E R is the seismic energy, T 0 the time window of 1 min, V 0 the average S-wave velocity of 3.5 km/s, r i the hypocentral distance, ρ the average density of 2700 kg/ m 3 , A i the square root of the sum of squared velocity amplitude of the three components in the 2-8 Hz band, t i the travel time as calculated from the hypocenter and the average S-wave velocity, f c the center frequency of 5 Hz, and Q the quality factor. We adopted Q values, which depend on the hypocentral distance, reported by Yabe and Ide (2014) and corrected for the effect of the free surface.

Detectability of tectonic tremor and redefinition of ETS events
In this study, we typically detected about three times as many tremor hypocenters as listed in the NIED hybrid catalog, making it possible to analyze the detailed growth process of an ETS event. We located 52,776 tremor hypocenters during 20 ETS events from our semblance analysis of the array data. The spatio-temporal distribution of tremors is in overall agreement with the result of Sagae et al. (2021), which was based on the MUSIC high-resolution f-k method, although their result had a time resolution as fine as 10 s. Figure 3 shows the space-time plot of tremor in ETS event number 17 (ETS 17), during 9-14 January 2014. The temporal variation in tremor activity we observed is consistent with that of the NIED catalog. The large number of hypocenters enabled us to recognize bilateral migration of the hypocenters along strike during ETS 17 (Fig. 3). The migration velocity of the main front of tremors is approximately 10 km/day, consistent with the typical case in the Kii Peninsula (Obara 2010). It is usually difficult to distinguish whether tremors belong to an ETS event or are background activity, causing uncertainty in the spatio-temporal range of an ETS event. To examine the growth process of an ETS event, we redefined ETS events on the basis of hypocenter clusters. We defined clusters as sets of tremor events that occurred within 1 h and 10 km epicentral distance, and the cluster hypocenter was defined as the centroid of a cluster comprising more than 10 tremors. The start and end points of an ETS event were defined as the centroids of the first and last clusters, and these points in turn defined the duration of an ETS event (Table 2). Three ETS events (ETSs 1, 8, and 11) contained too few clusters with centroids for accurate estimates, and were excluded from the following analyses and discussion.
As mentioned previously, short-term SSEs were geodetically detected in the analyzed period of this study. Table 3 lists the parameters and Additional file 2: Fig.  S2 shows the locations of faults associated with these short-term SSEs that were related to the ETS events of   Table 3 Short-term slow slip events estimated from geodetic data and their corresponding ETS numbers Geodetic data consisted of strainmeter, tiltmeter, and groundwater pressure records collected by GSJ and NIED. Seismic moments were calculated using a rigidity of 40 GPa. Asterisks denote SSEs outside the ETS periods determined in this study (Table 2); most SSEs are included within the semblance analysis periods (Table 1). Fault models for ETS 2 and ETS 3 were estimated independently by GSJ and NIED Date (yyyymmddhh) Duration  (Table 2). The fault parameters for ETSs 2 and 3 have been estimated independently by GSJ ) and NIED (Hirose and Kimura 2012). Despite the fact that estimations were constrained to rectangular faults and uniform slip, for most ETS events the locations of these faults appeared to be consistent with the distribution and duration of tremor. Geodetic estimations for ETSs 7, 8, 13, 14, and 20 indicated inconsistent duration of short-term SSEs with that of ETSs (Tables 2 and  3), although their durations were approximately defined by the occurrence of tremor within their respective semblance analysis periods (Table 1), except for the third events of ETS 7. Therefore, the uncertainty may be a consequence of how the start and end times of an ETS event were defined. For ETS 7, the eastern half of the fault of the second event and the entire fault of the third event were outside the area where tremor was detectable. For the sixth event of ETS 20, GSJ detected tremors during this event to the southwest of the fault of the sixth event , outside the area where we could detect tremor. It is noteworthy that during ETS 20, the locations of tremor and fault activity appeared to migrate in similar trajectories ( Fig. 5 and Additional file 2: Fig. S2). The ratio between the total radiated energy (Table 2) and the geodetically estimated seismic moment (Table 3) for the ETS events ranged from 10 -11 to 10 -9 , consistent with similar determinations for very-low-frequency earthquakes ) and SSEs (Maeda and Obara 2009).

Stationary patches of high seismic energy release
In this study, we used the radiated energy of tremor in the 2-8 Hz band as a proxy of the seismic moment of slow slip in a source region; likewise we used the radiated energy rate as a proxy of the seismic moment rate. Figure 4 plots our estimates of the cumulative energy of tremor during ETS 3 along with tilt records at Hi-net station KRTH during ETS 3 as reported by Hirose and Kimura (2012). Tidal effects and atmospheric pressure response were corrected with BAYTAP-G (Tamura et al. 1991) and detrended. The slope of the cumulative energy is gentle from 12 to 13 September, steep from 14 to the middle of 15 September, and gentle afterward. The same pattern is evident in the tilt records. This consistency between the two observations suggests that the energy released during tremor events may be associated with the slip rate during SSEs (Maeda and Obara 2009), and that the radiated energy rate of tremor can be used as a proxy of the moment rate of slow slip Rubin and Armbruster 2013). Kano et al. (2018) reported that the slip rate of slow slip and the radiated energy rate of tremor vary in tandem in western Shikoku. If this relationship also exists in the entire Nankai subduction zone, then the areas of high tremor energy release found in this study could be interpreted as areas of high slip during SSEs, i.e., areas of high seismic moment release. The consistency of the large temporal changes in both cumulative energy and tilt records in this study also suggests that most of the slow-slip seismic moment release occurred in areas of relatively high energy release. The strong spatial correlation between typical energy release rates of tremor and cumulative slip of SSEs in the Kii Peninsula (Yabe and Ide 2014) supports our use of radiated energy as a proxy of the seismic moment of slow slip in the study region.
On the other hand, previous works have noted that slip may occur without tremor during a short-term SSE (Obara et al. 2011;Wech and Bartlow 2014;Michel et al. 2019a). If tremorless slip had been dominant in the ETS events we analyzed, then energy release by tremor could not be a proxy for slow slip during the ETS events. However, the close correspondence of large tilt changes with energy release by tremor, as shown for example in Fig. 4, implies that tremorless slip neither dominated the seismic moment release by slow slip nor strongly influenced the growth process of most of the ETS events in this study. Figure 5 depicts the migration of tremor hypocenter clusters and the distribution of the cumulative seismic energy of tremors for three ETS events with different growth patterns (the other 14 ETS events are shown in  (Table 3) Additional file 3: Fig. S3). ETS 3 displayed unilateral migration from the southern end, ETS 20 displayed unilateral migration from the northern end, and ETS 9 displayed bilateral migration starting from the central part (upper panels in Fig. 5). All three events displayed patches with relatively high cumulative seismic energy release ('high-energy patches' hereafter) that appeared consistently in the same locations, regardless of the starting location and the migration direction (lower panels in Fig. 5). Therefore, the location of these patches appears to be controlled by physical or frictional properties of the fault plane that are inherent to the plate interface. The existence of asperities that rupture repeatedly during multiple ETS events is consistent with the presence of such asperities estimated from geodetic inversion in Shikoku  and documented in LFE epicenters in northern Cascadia .
Velocity structure heterogeneities may cause tremor to concentrate in certain locations. Among hypocenters determined by seismic tomography by Matsubara et al. (2019), those determined within a one-dimensional velocity structure differ from those determined within a three-dimensional structure by 1 km on average, and by as much as 8 km for hypocenters at depths of 0-50 km in this region (M. Matsubara, personal communication). This difference among hypocenters is small relative to the size of the high-energy patches. The effect of a heterogeneous velocity structure on the tremor locations in this study would be comparable to that of Matsubara et al. (2019). It is therefore unlikely that a heterogeneous velocity structure accounts for the observed high-energy patches. Ghosh et al. (2009) reported that the cumulative moment of tremor in an ETS event in Cascadia was spatially heterogeneous and characterized by several patches of high moment release. The distribution of these patches reported by Ghosh et al. (2009) appeared to be complementary to the slip surfaces of regular earthquakes, suggesting that inherent properties on the plate interface control the patch distribution. Figure 6 shows relationships between the location of high-energy patches, start points, and along-strike rupture extent (defined by tremor cluster hypocenters) of the 17 ETS events (see also Additional file 3: Fig. S4). It is notable that nearly all of the start points are outside the high-energy patches. This feature is similar to  et al. Earth, Planets and Space (2021) 73:59 the relationship between hypocenters and asperities of regular earthquakes (e.g., Yamanaka and Kikuchi 2004;Mai et al. 2005) and suggests that the rupture of an ETS event initiates in a low-strength area. This implies that the higher-energy patches are areas with relatively high frictional strength. It is also notable that start points tend to be on the deeper, down-dip side of the plate interface, whereas the high-energy patches are on the up-dip side, implying a depth-dependent variation in frictional strength on the plate interface Sweet et al. 2019). The down-dip concentration of the start point is consistent with the observations of major tremors in the Kii Peninsula reported by Yabe and Ide (2014). Moreover, the along-strike distribution of start points appears to be heterogeneous. Some start points are concentrated on the northeastern side of patch A. This concentration, together with the rupture of patches A and B, causes the dominantly southwestward migration of tremor in the analyzed area (Obara 2010). Focusing on the rupture extent of ETS events, we found that ETS events tended to terminate in or around highenergy patches. This suggests that patches have relatively high frictional strength. In that case, the complete rupture of a high-energy patch would extend the growth of an ETS event. On the other hand, the termination of the rupture in or around a patch would prevent the further growth of an ETS event.

Relationship between start point, rupture extent, and high-energy patches
We have mentioned the role of high-energy patches, or tremor asperities (Ghosh et al. 2012), in the rupture propagation of ETS events. It may appear unlikely that tremor asperities actively influence the growth of ETS events because the total seismic moment of tremors and LFEs is orders of magnitude smaller than the seismic moment of a contemporaneous SSE (Kao et al. 2010), and tremor and LFEs are generally interpreted as a passive manifestation of underlying slow slip. However, some studies have argued for an active role of tremor or LFEs in the propagation of slow slip. Shelly (2015) found that on the deep San Andreas fault, tremor propagates most effectively through regions of abundant tremor production, and does not propagate through large gaps in tremor production, implying an active role for tremor in slip propagation under the interpretation of rapid tremor migration as a self-regulating cascade of seismic ruptures along the fault. Rubin and Armbruster (2013) suggested that tremor could influence the evolution of slow slip at small scale lengths in their analysis of radiated energy rates and tremor migration in a secondary slow slip front in a 10 km 2 area of concentrated tremor in Cascadia. These studies restricted the active role of tremor to the small-scale propagation of slow slip.
However, tremor may play an active role not only in small-scale slow slip propagation, but also during entire ETS events. Yabe and Ide (2014) compared the migration pattern and seismic energy release rate of tremor, which they referred to as the strength of tremor patches, in the Nankai, Cascadia, Jalisco, and South Chile subduction zones and suggested that tremor does control the occurrence of SSEs to some extent. They proposed a qualitative model, an expansion of the model of Ando et al. (2012) which used the spatial distribution of variations in the strength of tremor patches to explain the behavior of slow earthquakes and tremor. For regions of heterogeneous tremor-patch strength, such as the Kii Peninsula, Yabe and Ide (2014) suggested that the rupture of weak tremor patches transfers the stress accumulated by tectonic loading to strong tremor patches, and that strong tremor patches serve as switches that initiate large SSEs. This model explains the observed characteristics of the ETS events in this study: their initiation outside high-energy patches, their termination in and near high-energy patches, and their cascading growth discussed below. The role of tremor asperities in the growth of SSEs was also examined by the simulation studies of Ampuero (2014, 2017), who developed mechanical models consisting of frictionally unstable asperities embedded in a frictionally stable fault zone matrix. Luo and Ampuero (2017) used numerical simulations to investigate two end-member models, the SSE-driven-tremor model and the tremordriven-SSE model, and concluded that although both models quantitatively reproduce the observed characteristics, the latter model is more plausible because it is based on a more conventional friction law and better reproduces the observed characteristics of SSE and tremor without fine-tuning of the model parameters. These studies suggest that high-energy patches can play an active role in the growth of the ETS events reported in this study. Figure 7 shows two examples of rapid tremor reversal (RTR) (Houston et al. 2011) observed in this study. During ETS 2, the main front of tremor propagated from patch A to patch B and then propagated back to patch A. During ETS 20, the main front of tremor propagated Fig. 7 Evolution of RTRs during a ETS 2 and b ETS 20. The top three panels show the cumulative seismic energy release during the period of each snapshot (upper left corner). Black rectangles show the locations of high-energy patches on the subducting Philippine Sea plate interface. The second plot from the bottom shows the spatio-temporal distribution of tremor hypocenters, with colors corresponding to radiated energy rate and dotted lines indicating migration speeds. The bottom plot shows the tidal shear stress at a depth of 35 km on a fault plane (with strike, dip, and rake taken from corresponding short-term SSEs listed in Table 3) below patch B for ETS 2 and below patch C for ETS 20. Positive shear stress promotes slip on a fault from patch B to patch C, then back to patch B. It appears in both cases that a high-energy patch was reruptured when RTR occurred. Our observations show clearly that the patches release a large portion of their seismic energy in the first rupture and a smaller portion in the second rupture, accompanied by RTRs (Fig. 7). The migration speed during RTR reached approximately 250 km/day, which is of the same order as RTR velocities observed previously in the Kii Peninsula ) and in Cascadia (Houston et al. 2011), whereas the migration speed during the main rupture front (i.e., the first rupture of the patches) was approximately 15 km/day, similar to speeds of the main rupture front in other ETS events in this region. This migration speed is also comparable to that estimated by Obara (2010); 2-20 km/day and averaging 10 km/day.

Relationship between rapid tremor reversal and high-energy patches
The relatively high energy release rates of tremor during RTRs observed in this study are consistent with previously reported high tremor amplitudes during RTRs (Thomas et al. 2013;Rubin and Armbruster 2013;Bostock et al. 2015;Peng et al. 2015). The total radiated energy of tremor during the RTR is estimated to be roughly 7 × 10 6 J for ETS 2 and 6 × 10 6 J for ETS 20, or 4% and 2% of the entire radiated energy, respectively. This small fraction of radiated energy related to RTR implies that the fraction of the total seismic moment of slow slip related to RTR is similarly small, as observed by Hawthorne et al. (2016), if the same proportionality holds between the radiated energy rate of tremor and the seismic moment rate of slow slip.
Large-amplitude RTRs are often tidally modulated Thomas et al. 2013;Royer et al. 2015), whereas others are triggered by a cascading failure of brittle asperities (Peng and Rubin 2016). To investigate the cause of the RTRs in ETS 2 and ETS 20, we used the TidalStrain.2 code (Hirose et al. 2019) to calculate the tidal stresses at a depth of 35 km below patch B for ETS 2 and below patch C for ETS 20 because the RTRs appear to initiate from these patches. The assumed fault geometry and slip (strike, dip, and rake) were taken from the parameters of a corresponding SSE (Table 3). We adopted shear stress (∆τ) as the index of tidal stress because the frictional coefficient is close to 0 for tremor (Houston 2015). The timing of the RTRs does not coincide with a peak in shear stress or a positive shear stressing rate ( Fig. 7), ruling out modulation of these RTRs by tidal shear stress as reported by Rubin and Armbruster (2013). This result agrees with that of Yabe et al. (2015), who ruled out such a tidal modulation in the Nankai subduction zones. Instead, these RTRs appear to have been triggered by the failure of high-energy patches in the absence of high-energy rate tremor. The trigger of the RTRs without high-energy rate tremor is consistent with the observation of RTRs of Peng and Rubin (2016) and the simulation of Ampuero (2014, 2017).

Cascading growth process of ETS events
Comparing the temporal evolution of ETS events with the same start points effectively removes one possible confounding variable in the growth process of ETS events. Figure 8 displays the temporal evolution of three ETS events (ETS 9, 12 and 17) that initiated only a few kilometers apart. These events started at the down-dip end, and the rupture propagated up dip at about 1 km/h for about 12 h, which is consistent with the migration speed of 1 km/h reported for these events by Sagae et al. (2021). ETS 12 ruptured patch B and terminated without rupturing patches A and C, whereas ETSs 9 and 17 ruptured patch B and then propagated bilaterally to patches A and C. In ETS 17, patch A ruptured before patch C. The rupture of ETS 17 continued longer than that of ETS 9, and the radiated seismic energy from patches A and C was larger in ETS 17 (Fig. 8b). These results suggest that differences in ETS growth are controlled by the rupture of strong patches, as suggested by Yabe and Ide (2014), even for ETS events with nearly identical start points.
For regular earthquakes, the existence of a nucleation phase, or slow initial phase, has been controversial, and some studies have reported that the eventual size of an earthquake is scaled as the cube of the duration of the nucleation phase (e.g., Ellsworth and Beroza 1995). We investigated whether such a phase can be found for ETS events, again using the trio of ETS 9, 12, and 17 (Fig. 8c). Given the possible uncertainty of the start time of ETS events, we assigned a consistent start time for the three events by aligning their energy release curves at the beginning of the rupture of patch B, then considered the relation between energy release during the day before and the day after this start time (Fig. 8c). We also assigned the start time as the time of the first (See figure on next page.) Fig. 8 a Evolution of the energy release for three ETS events with nearly identical start points. Each snapshot encompasses the time period shown in the upper left corner and displays the cumulative seismic energy release during that period. Black rectangles represent the locations of high-energy patches on the subducting Philippine Sea plate interface. Red stars indicate the start point of the ETS events. b Cumulative seismic energy release curves of the three ETS events. c, d Comparison of cumulative seismic energy release during the initial parts of the three ETS events in which c the start of the rupture in patch B and d the time of the first tremor event within the period of semblance analysis (Table 1) serves as the zero point on which the curves are aligned so as to better compare the earliest phase of the ETS events Nakamoto et al. Earth, Planets and Space (2021) 73:59 tremor event within the period of the semblance analysis (Table 1) under the assumption that all tremor events we detected were related to ETS events (Fig. 8d). Despite the large difference in their sizes (Fig. 8b), the three events had very similar initial phases before the start time; that is, there was no size-dependent growth process evident in the initial part of the ETS events. This result suggests that the growth of ETS events is controlled by how highenergy patches rupture, and suggests that ETS events have a cascading rupture process similar to that reported for regular earthquakes (Ide 2019). A logarithmic plot of the data shown in Fig. 8b (Fig. 9a) allows us to focus on the long-term cumulative energy curve of ETS events after 10 4 s (about 2.8 h), a time range that emphasizes variations between events. The cumulative energy of slow earthquakes is proposed to be proportional to the seismic moment with a constant of approximately 10 -10 , and the cumulative energy curves in Fig. 9a appear roughly consistent with a scaling relationship between seismic moment (M O ) and duration (T) of M O ∝ T for slow earthquakes (Ide et al. 2007), shown as the pink stripe in Fig. 9. However, it is apparent that the slopes of the cumulative energy curves are not constant. The slope is steep for the time interval 2-4 × 10 4 s to 5-9 × 10 4 s, when the cumulative energy (E) scales approximately as the cube of T, or E ∝ T 3 . The slope becomes gentler after 5-9 × 10 4 s, when the cumulative energy scales approximately as the lapse time, E ∝ T (Fig. 9a). The first time interval corresponds to the rupture of patch B, and the second corresponds to the rupture propagation along strike from patch B to patch A or C. Therefore, we suggest that ETS events might grow as M O ∝ T 3 during rupture expansion in a high-energy patch and as M O ∝ T during rupture propagation along strike within the ETS zone. This transition of the scaling relation due to the bounded or one-dimensional growth of SSEs within the ETS zone was also suggested by Gomberg et al. (2016). In this case, the scaling relation of M O ∝ T 3 for SSEs in Cascadia reported by Michel et al. (2019b) might indicate the unbounded growth of a SSE within the ETS zone. Moreover, the transition we found in the scaling relation between M O and T differs from the result of a numerical simulation of slow earthquakes based on a Brownian walk model, which shows a transition from M O ∝ T 2 to M O ∝ T (Ide 2008), although this difference might reflect fluctuations in the temporal evolution of ETS events. The cumulative energy curves for all 17 ETS events in this study are quite scattered (Fig. 9b), and they appear to have no common evolution process signified by a systematic transition of the scaling between E and T. On the other hand, most of the events approximate the M O ∝ T scaling relationship (Ide et al. 2007) when considered as whole events (Fig. 9b). This could be consistent with the M O ∝ T 0.811 scaling relationship for 61 SSEs in Shikoku in the Nankai subduction zone (Hirose and Kimura 2020).

Conclusions
We analyzed seismic array records of tectonic tremor beneath the Kii Peninsula in the Nankai subduction zone to trace the detailed growth of ETS events. Our examination of tremor in 17 ETS events revealed the existence of stationary patches on the plate interface that release high radiated seismic energy. The location of these Fig. 9 Logarithmic plot of the evolution of cumulative seismic energy release (a) for the three ETS events shown in Fig. 8b and b for all ETS events. The pink zone represents a scaling relation of M O ∝ T (Ide et al. 2007) with a proportional relation to seismic energy with a constant of 10 -10 . Black lines indicate scaling relations of E ∝ T, E ∝ T 2 , and E ∝ T 3 high-energy patches is independent of the start point and migration pattern of ETS events, suggesting that the distribution of the patches reflects heterogeneities in frictional properties that are inherent to the plate interface. The start points generally lie outside the high-energy patches, similar to the spatial relationship between hypocenters and asperities of regular earthquakes. In some ETS events, the high-energy patches have a role in terminating the rupture; in others, the high-energy patches rupture twice in an RTR-like pattern. We also compared the initial part of ETS events with different sizes but nearly identical start points, finding that the final size of ETS events appears to be controlled by the extent of the rupture of the high-energy patches; that is, the final size of an ETS event cannot be estimated from its initial growth phase. These features indicate that ETS events have an essentially cascading rupture process similar to that of regular earthquakes. The growth of ETS events appears to be controlled by how they rupture highenergy patches.