Aseismic slip and seismic swarms leading up to the 2024 M7.3 Hualien earthquake

Domain: earth_science
Published: 2025-10-13
Editorial outcome: major_revision

Abstract
Understanding the role of aseismic slip in earthquake cycles is essential for assessing seismic hazards and short-term forecasting. Eastern Taiwan’s double-vergence suture zone, where the Philippine Sea Plate subducts beneath the Eurasian Plate, experiences frequent M≥6 earthquakes and widespread aseismic slip, making it an ideal natural setting to study earthquake triggering processes. Here we demonstrate how aseismic deformation contributed to the April 3, 2024 Mw 7.3 Hualien earthquake by analyzing a 24-year catalog of repeating earthquake sequences (RESs) and earthquake swarms. We find that six out of nine swarms in the epicentral area, northern Longitudinal Valley, were accompanied by increasing aseismic slip rates, as revealed by RESs on the west-dipping Central Range Fault (CRF). A notable aseismic slip episode in 2021 indicated by GNSS signals, the accelerated RESs-derived slip rate, and a four-month-long swarm sequence with high diffusivity (~5.2 m²/s), suggests joint contributions from over-pressured fluids and deep fault creep. Following this episode, a sequence of M6+ events occurred in 2022, and both seismicity and aseismic slip gradually increased again starting in 2023. Coulomb stress modeling indicates that cumulative aseismic and seismic slips since 2021 generated up to ~30 kPa positive stress on the eventual 2024 rupture, promoting fault weakening and shallower seismicity. This study provides compelling evidence for aseismic-slip-induced stress triggering of a major earthquake and highlights the importance of integrating aseismic processes into earthquake hazard models for collisional fault systems.

Introduction
Understanding the interaction between slow slip and regular earthquakes is key to advancing our knowledge of earthquake nucleation processes e.g., refs. 1,2. While some studies suggest that slow slip events (SSEs) may trigger large earthquakes by modifying local stress conditions e.g., refs. 3–6, others propose that SSEs relieve accumulated stress and reduce the likelihood of large events 7,8. Simulations of SSEs 9 indicate that recurrence intervals of long- and short-term SSEs may shorten prior to a major earthquake. However, a global catalog comparison of SSEs and nearby seismicity does not reveal a consistent correlation between slow slip and subsequent earthquake activity 10. These contrasting interpretations may arise from limited observation period, low detection resolution, or the inherently variable behavior of slow slip phenomena. To date, compelling evidence of ubiquitous precursory slow slip remains limited.

In this study, we examine aseismic slip behavior in eastern Taiwan using two complementary indicators: repeating earthquake sequences (RESs) and earthquake swarms. Repeating earthquakes, which reflect localized fault creep, are sensitive to variations in aseismic loading e.g., ref. 6. Earthquake swarms, meanwhile, often indicate transient aseismic slip and/or fluid overpressure, and are highly sensitive to subtle changes in stress and pore pressure e.g., refs. 11–13. Laboratory and theoretical studies show that changes in pore pressure can reduce effective normal stress, promoting a spectrum of slip modes including slow and fast earthquakes 14–16. Taiwan provides an exceptional natural laboratory to investigate these processes due to its rapid tectonic deformation, frequent large earthquakes, and dense seismic instrumentation. Here, we explore the interaction between aseismic and seismic slip preceding the April 3, 2024 Mw 7.3 Hualien earthquake.

Taiwan is one of the most seismically active regions globally, especially along the active Longitudinal Valley (LV) suture zone in eastern Taiwan. This zone, resulting from the collision between the Luzon Arc on the Philippine Sea Plate and the continental crust of the Eurasian Plate, hosts two parallel, head-to-head fault structures with opposing dips: the east-dipping Longitudinal Valley Fault (LVF) to the east and the west-dipping Central Range Fault (CRF) to the west e.g., refs. 17–19. Since 2000, these fault systems have produced 13 M≥6 events (Table S1), revealing complex interaction patterns, including (1) simultaneous ruptures or delayed triggering during sequences of M7 events e.g., refs. 20,21 and (2) out-of-phase occurrences of M≥6 events due to the stress shadow effect 22. Among these, two events exceeded magnitude 7.0, as shown by numbers 10 and 11 in Fig. 1 and Table S1 of the supplementary material.

The April 3, 2024 Mw 7.3 Hualien earthquake likely nucleated on the LVF 23, as inferred from the alignment of the earlier stage of aftershocks (see Movie S1 of the supplementary material). This event was followed by two M6-7 and 25 M5-6 aftershocks within a month, which reveal conjugate ruptures on both the northern LVF and CRF 24 (purple circles in Fig. 1b). In contrast, the September 18, 2022 Mw 7.0 Guanshan-Chihshang earthquake sequence initiated on the southern CRF and triggered aftershocks on both the CRF and the creeping segment of LVF, the Chihshang Fault (marked by light purple circles in Fig. 1b,d) 21,22,25. The seismicity, aftershocks, and slow slips along both creeping and locked segments of the LVF and CRF present a unique opportunity to improve understanding of fault interactions via aseismic slip and the role of aseismic slip in earthquake cycles.

In this study, we compile a 24-yr catalog of RESs and earthquake swarms to assess the spatiotemporal association between slow slip indicators and regional seismicity, and further, investigate aseismic slip evolution leading up to the 2024 Mw 7.3 Hualien mainshock. In the following sections, we (1) quantify the aseismic slip rate variations using RESs data, (2) analyze swarm migration and dynamics, (3) identify GNSS-observed slow slip episodes, and (4) evaluate seismicity rate changes preceding the mainshock. Finally, we use Coulomb stress modeling to explore how these aseismic processes may have contributed to triggering the 2024 Mw 7.3 event.

Results
Aseismic slip rates from repeating earthquakes
RESs involve groups of events characterized by nearly identical waveforms, magnitudes, and locations, representing repeated ruptures of the same fault patches. Their short recurrence intervals imply rapid loading from surrounding aseismic slip, making them effective indicators of fault creep and interseismic deformation rates 26–29. In eastern Taiwan, RESs have been documented along the creeping segments of the LVF and CRF 30–34. We updated the previously established RESs catalog 30 from 2000–2011 to 2012–2024, to include events through May 2024. The RES detection approach is described in Methods.

The updated RESs catalog contains 649 events across 148 sequences from January 1, 2012–April 18, 2024, with magnitudes ranging from 2.0 to 5.1. Merging this with the January 1, 2000 through December 31, 2011 RES catalog 30, the integrated dataset comprises 1499 RES events in the 148 sequences. Each sequence includes 3–34 events. As shown by diamonds in Fig. 1, these events cluster into two spatial groups: in the northern LV (north of 23.6°N), RESs occur on the west-dipping CRF, while in the southern LV (south of 23.6°N), they cluster along the east-dipping LVF.

Using the scaling relationships developed for creeping sections in eastern Taiwan 30, the slip for each repeating event can be obtained using the scaling relationship between the slip estimate (d) and seismic moment (Mo) as log d = −1.21 + 0.11 log Mo and log d = −1.96 + 0.14 log Mo for the LVF and CRF, respectively (see Methods for how these scaling laws were inferred). Slip histories from multiple RESs were then combined to estimate the spatiotemporal variations in aseismic slip driving the repeating events. The regional cumulative slip was calculated in one-day increments by measuring the change in slip between consecutive days. To ensure that the calculation window exceeded the length of any data gaps but remained shorter than the targeted recurrence pattern, we applied a 180-day averaging window. This smoothing reduces short-term variability and enables a more robust determination of the slip-rate time series 30.

Figure 2 presents the temporal distribution of RESs events and the corresponding aseismic slip rate, while Fig. S1 in the supplementary material visualizes the along-strike spatiotemporal evolution of M≥2 seismicity density and RESs-derived aseismic slip. We systematically tested the sensitivity of the computed slip rate to different averaging windows (Fig. S1b–c). Shorter averaging windows yield higher temporal resolution but are prone to noise when below the RESs inter-event time (mean: 15 days; maximum: 127 days). Longer windows increase reliability but may introduce time lags 30. Near-zero slip rates using 1-month time window (Fig. S1b–c) typically reflect gaps in RESs detection.

In the southern LV area where RESs are widely distributed along the east-dipping LVF, from near surface to 30 km (Fig. 1d), the highest aseismic slip rate peak occurred shortly after the 2003 Mw 6.8 Chihshang earthquake, coinciding with a sharp increase in RES activity. Following the Mw 6.8 mainshock, the RESs occurrence rate tripled from 13.7 to 41.6 events per year, which persisted for 6 years. In the northern LV where the RESs mainly occurred on the west-dipping CRF below the depth of 10 km, the RESs occurrence rate is generally higher than in the southern LV; The small fluctuations in aseismic slip rate tend to show annual variation 30. The most pronounced CRF aseismic slip acceleration followed the April 18, 2019 Mw 6.1 event, with RESs activity increasing from 33.3 to 45.8 events/year.

Not all large earthquakes in eastern Taiwan influence aseismic slip as captured by RESs. Along the northern CRF, where more than five M≥6 earthquakes occurred during the study period, the RESs-inferred aseismic slip rate remains relatively stable. This may reflect the fact that RESs tend to occur at depths of 10–25 km, whereas most large earthquakes ruptured shallower parts of the fault (<10 km). Only RESs located in close proximity to mainshock hypocenters (e.g., Number 6 in Fig. 1) show evidence of stress-driven acceleration in aseismic slip.

Aseismic slip behavior from earthquake swarms
Aseismic slip can also manifest through earthquake swarms, clusters of earthquakes lacking a clear mainshock and atypical aftershock decay pattern. These swarms are generally associated with fluid intrusions and elevated pore pressure e.g., refs. 35–37, which can induce fault unclamping and promote slow slip 38,39. A previous study in Taiwan 26 analyzed 153 M≥3 swarm sequences (4726 events) and their relationship with 59 M≥6 earthquakes, finding limited evidence of precursory swarms. Instead, swarms more often followed large earthquakes. As demonstrated by the previous study 26, a significantly higher proportion of swarm events followed M≥6 earthquakes when the separation distance was less than 10 km, highlighting a spatial dependence in post-mainshock swarm triggering.

To build a more comprehensive swarm catalog, we applied a composite clustering method 26 that integrates multiple declustering algorithms. Swarm sequences were identified from events with M≥2.0, consistent with the estimated magnitude of completeness (Mc) for this region (see Methods for details). From January 1, 2000 to April 30, 2024, we identified 15 swarm sequences comprising 4109 events with magnitude ranging from M 2.0 to 5.9 (Table S2 and Fig. S2). Of these, 82% of the swarm events (in 10 sequences) occurred in the northern LV (Fig. S2a), coinciding with the aftershock zone of the 2024 Mw 7.3 Hualien earthquake and along the west-dipping CRF (Fig. 1b,c, Fig. S2c). These northern swarms are spatially associated with high Vp/Vs zones (cross-section in Fig. S2b–c), suggesting a fluid-saturated environment. In addition, the composite method yields a declustered earthquake catalog that excludes mainshock-aftershock sequences and swarm-like clusters, which is used for background seismicity analyses.

The longest swarm sequence (Sequence 9) occurred from April 18 to August 18, 2021, comprising 796 events, including nine 6>M≥5 earthquakes (Fig. 3). Approximately 2 months after the first event, 86% of swarm events occurred at depths shallower than 15 km, ultimately reaching regions of high Vp/Vs ratio and low Vp and Vs 27. The swarms mainly spanned depths of 5–22 km, partially overlapping with RESs depths (13–22 km, green squares in Fig. S11). This swarm sequence demonstrated bilateral migration, with 69% of events propagating upward and 31% downward. Compared to the vertical direction, northward migration was less pronounced (Fig. 3a, b and Movie S2).

To quantify the swarm’s migration, we applied a diffusion model to the swarm triggering front. Assuming constant permeability and pore pressure conditions, the hydraulic diffusivity (D) was estimated by fitting the diffusion equation r = \sqrt{4\pi Dt}, where r is the distance from the first event in a swarm sequence to the triggering front, t is the elapsed time from the start of diffusion, and D is the hydraulic diffusivity 28. The best-fit curve was selected to span the majority of 90th percentile distance points (pink squares in Fig. 3e) obtained every 20 earthquakes 29. The minimum RMS misfit between the observed and predicted triggering front is obtained using moving time bins containing 20 events, with 10 events overlapping between consecutive bins.

A sensitivity test was performed by sequentially removing early data points (the first 1 to 20 pink squares in Fig. 3e). In other words, the first time we removed the first point (90th percentile obtained using the first 20 events) to measure D1, the second time we removed the initial two points for D2, and the third time D3 for removing three points, and finally, D20 for removing the first 20 points at the earlier stage of the sequence. This yielded a distribution of diffusivity estimates (D1 to D20), from which we computed the mean and standard deviation. For Seq. 9, the mean hydraulic diffusivity was approximately 5.4 m²/s, with standard deviation of 1.1 m²/s, indicating rapid fluid migration from the mid-crust.

Key parameters for each swarm sequence in northern LV area are listed in Table S2, including onset time, duration, sequence duration, number of M≥5 events, seismic moment (maximum and cumulative), b-value, mean diffusivity (D), and standard deviation (SD). Notably, a new swarm initiated 19 days after the 2024 Mw 7.3 mainshock on April 22, was excluded from this study due to the catalog cutoff at April 30, 2024. We found that across all swarm sequences, the cumulative moment scales positively with both maximum event size and swarm duration, but inversely with b-value. All sequences except Seq. 1 had b-values < 0.91. This pattern can be explained by the tendency for sequences with a larger maximum event to also include a higher proportion of relatively large events, resulting in lower b-values. Such sequences, often led by relatively large events (e.g., M>5), tend to persist for longer durations and release greater cumulative seismic moment.

Migration was observed in all swarm sequences but showed no consistent directivity (Figs. S3–S11). High D values were generally associated with greater uncertainty, suggesting that diffusivity estimates exceeding ~10 m²/s should be treated with caution. In contrast, focused swarms (e.g., Seqs. 8 and 9) featured more events and yielded lower diffusivity standard deviations, providing more reliable estimates.

Spatiotemporal association between swarm versus repeating earthquakes
We examined the temporal relationship between swarm sequences and deeper aseismic slip using the RESs time series. Figure 4c compares individual swarm timelines with the evolving aseismic slip rate. To assess potential triggering or causal relationship, we aligned all swarms by their onset (time = 0) and stacked the surrounding aseismic slip rate histories across a ±180-day window (Fig. 4b). We found that 6 out of 9 swarm sequences were preceded by an increase in aseismic slip rate, suggesting a possible causal link between deep aseismic slip and the initiation of swarm activity. Exceptions include Seqs. 4-6, where no clear precursory acceleration was observed.

To test the robustness of the temporal relationship between the changes in aseismic slip rate and swarm activity, we performed a sensitivity analysis by calculating the average aseismic slip rate before (dpre) and after (dpost) each swarm onset across sliding time windows of 1 to 12 months. The difference (df = dpost - dpre) quantifies whether slip rate accelerated before or after the swarm onset. A negative df indicates precursory acceleration, while a positive df suggests slip increase following the swarm. Despite some variability across time windows, Seqs. 3, 4, 5, and 6 consistently exhibited negative df values, while Seqs. 2, 7, and 9 showed positive df (Fig. S12). Most consistent results were obtained using 3- to 6-month averaging windows, which revealed precursory acceleration in 6 of 9 cases, supporting the hypothesis that increased aseismic slip often precedes swarm activity.

Seqs. 2 and 9 show the strongest temporal and spatial alignment between swarms and RESs (i.e., relatively large number of repeating events are activated in the close neighborhood of swarm events, indicated by green squares vs. colored circles in Figs. S4 and S11), with significant slip rate changes near swarm onsets (Fig. S12).

To further examine the spatiotemporal relationship between RESs accelerations, swarms, and M≥6 events, we plotted separation distance against time difference between each pair of phenomena in Fig. S13. Here, the aseismic slip acceleration is defined by the peak in RESs-derived slip rate history in Fig. 4c. (1) Swarm vs. RESs acceleration (Fig. S13a): Spatial separations range from 10 to 50 km, with time differences generally < 0.4 yr. However, smaller spatial separation does not consistently correspond to shorter interaction times, suggesting different underlying mechanisms. (2) Swarm vs. M≥6 events (Fig. S13b): A positive correlation is observed, as that closer swarm-mainshock pairs tend to occur with shorter separation times. The closest pairing (27 km) corresponds to a long-duration swarm and short time lag. (3) RESs acceleration vs. M≥6 events (Fig. S13c): No clear trend is observed, despite spatial separations of 5–40 km and time lags < 0.35 yr. These results indicate that while swarms and M≥6 events exhibit possible triggering relationships, RESs accelerations appear more loosely coupled with subsequent large earthquakes. Further statistical examination is needed to confirm and strengthen these correlations and interpretations.

Geodetic observation during the swarm sequences
To quantify crustal deformation associated with the swarm sequences, particularly after 2019 where RESs-derived slip rate began to accelerate (Fig. 2a), we analyzed GNSS data pre-processed by the Taiwan geodetic model at the Institute of Earth Science, Academia Sinica, Taiwan (https://tgm.earth.sinica.edu.tw/). To detect the transient deformation, we selected 11 GNSS stations (Fig. 5 and Fig. S14) and processed the time series data by correcting long-term linear trends, seasonal periodic motions, local M≥6 coseismic offsets and related post-seismic deformation (the GNSS processing methodology is detailed in the Methods section).

The corrected time series reveals transient deformation patterns coinciding with swarm activity (Fig. 5b, c). The 2019 swarm (denoted by red dots in Fig. 5a and lower panel of Fig. 5b, c), occurred approximately four months after the Mw 6.1 earthquake and ~30 km to its south, is found to produce ~10 mm of eastward and southward displacement at stations FLNM and DNFU. During the 2021 swarm sequence (denoted by black dots in Fig. 5a and lower panel of Fig. 5b, c) particularly between Mw 5.9 and Mw 5.5 events, gradual eastward displacements of ~15 mm were observed at stations YENL, SOFN, <10 mm at TUNM and HUAL, and up to ~20 mm at NDHU. We interpret these displacements as primarily reflecting a combination of post-seismic and aseismic deformation, rather than coseismic slip from the Mw 5.9 event at the depth of 15 km and its afterslip. This interpretation is supported by the fact that most seismicity between the Mw 5.9 and Mw 5.5 events was concentrated around 10 km depth (Fig. 3b, c), a depth unlikely to generate measurable surface displacements. Moreover, the GNSS displacement time series during this interval can be well modeled by a post-seismic exponential decay function without requiring a co-seismic step (Fig. S15).

To interpret the observed transient deformation, we performed a forward modeling assuming a west-dipping fault plane as the primary rupture surface. Due to the limited number of GNSS stations (only two or three) showing significant displacements near the swarm epicenters, the fault geometry was constrained simply using the locations and focal mechanisms of M≥5 swarm earthquakes. The dip and rake angles were fine-tuned to better fit the observations (see Table S3). Assuming that the observed transient deformation is solely tied to slip associated with the swarm events, we approximated the imposed fault slip based on the total seismic moment released during each sequence. For the 2019 swarm, this corresponds to ~2 cm of slip across a 15×20 km area. For the 2021 swarm (time period between Mw 5.9 and Mw 5.5 event), we consider two fault segments based on the depth distribution of swarm events: ~5 cm of slip on the shallow segment (40×10 km area) and ~4 cm on the deeper segment (40×18 km area), as illustrated by the blue vectors in Fig. S16.

The forward modeling results show that for the 2019 swarm, the predicted displacements account for only about half of the observed surface motion (red and black arrows in Fig. S16), suggesting a significant contribution from additional aseismic slip. In contrast, the predicted displacements for the 2021 swarm are slightly larger than the observations. These discrepancies may arise from simplifications in fault geometry and limited GNSS station coverage. A denser GNSS network and more detailed fault modeling will be necessary to better constrain the aseismic slip contributions.

Interplay of seismic and aseismic slip episodes in northern LV
In the region surrounding the active swarms (grey box in Fig. 4a), intense aseismic slip episodes were observed following the April 18, 2019, M6.1 earthquake. This is indicated by a sharp increase in the occurrence rate of RESs at the time of this M6.1 event (Fig. 2a), followed four months later by the onset of the 2019 swarm on August 2, 2019 (Seq. 8). Figure 6 summarizes the temporal evolution from 2019 to 2024, including swarm activity, RESs-derived aseismic slip rate, GNSS displacements corrected using monthly averages, and seismicity parameters.

In April 2019, an M6.1 earthquake occurred at the depth of 20 km in the northern part of the study area, spatially close to a cluster of deep RESs. These RESs responded with a notable acceleration in slip rate (first red dashed line in Fig. 6). Approximately 3.5 months later, a separate swarm sequence occurred to the south (first orange peak in Fig. 6a), accompanied by a discernible increase in the RESs activity, which is visible in the monthly average slip rate (green curve in Fig. 6b), along with ~10 mm of gradual displacement recorded at nearby GNSS stations (blue circles in Fig. 6e). Given the considerable spatial separation between the swarms, RESs, and the 2019 M6.1 event, we infer the presence of two distinct aseismic slip episodes: (1) a deep slip event at depths greater than 15 km in the north, where deep, widely distributed RESs aligned along the mountain strike responded to the nearby 2019 M6.1 event, leading to significant acceleration of aseismic slip rate; and (2) shallower slip episodes in the south, where aseismic slip was detected by nearby GNSS stations during the shallow swarm sequence.

In contrast, the prolonged 2021 swarm sequence from April 17 to August 29, spanned a broad depth range from ~24 km to near the surface and overlapped substantially with RESs (Fig. S11), yet occurred independently of any major earthquake. This episode marked the onset of the most consequential aseismic–seismic coupling observed in the 24-year record that may have initiated a sequence of fault processes that ultimately led to the 2024 Mw 7.3 Hualien earthquake. This swarm-RESs coupling, and its potential connection to the 2024 Mw 7.3 earthquake, requires detailed examination. The spatiotemporal evolution of this activity can be divided into four stages:

Stage 1 (April-August, 2021): Acceleration of aseismic slip and swarm activity. The intense earthquake swarm lasted from April 17 to August 29, initiated near 14 km depth and exhibited multiple waves of activity. During the first two months (first spike for 2021, Fig. 6a), swarm events clustered at depths of 15–20 km and migrated upward and northward. At the same time, the RESs-derived aseismic slip rate at 15–25 km depth began to accelerate, a trend most evident in the monthly average (green line in Fig. 6b). Increases in RESs’ moment rate and occurrence rate are also observed during this interval (Fig. 6b–d). In the later half of the sequence (second spike for 2021 in Fig. 6a), the swarm activity intensified at shallower depths, accompanied by a renewed rise in the RESs-derived slip rate. This indicates that deep aseismic deformation persisted while shallow brittle failure became more pronounced. Concurrently, GNSS stations near swarm zone recorded transient horizontal displacements (Fig. 6e).

Notably, the swarm period is marked by a sharp increase in shallow seismicity (≤15 km), reflected in both a greater number of events and elevated a- and b-values (Fig. 6f–h, dark blue). In contrast, the declustered background seismicity (Fig. 6i) shows no comparable change, indicating that the observed acceleration in shallow seismicity was driven primarily by the swarm itself rather than by background seismic processes. Over the full four-month period, the shallow seismicity rate increased from 4.39 to 10.75 events/day, while the deep seismicity rate remained steady at ~4 events/day (Fig. S17). The simultaneous acceleration of swarm activity, enhanced RESs-derived slip rate, and transient surface deformations suggest an aseismic slip episode that likely accompanied by fluid migration and activated both brittle and aseismic slip along multiple depth levels of CRF.

Stage 2 (September, 2021 to late-2022): Cascade-like M6+ events. After the 2021 swarm, both the RESs-derived slip rate and overall seismicity declined until late 2022. During this period, however, a sequence of large earthquakes occurred, including the March 22 M6.7 event in southern LV (number 7 in Fig. 1), the June 20 M6.1 event in northern LV (event 8 in Fig. 1), and September 17-18 M6.5 and M7.0 events in southern LV (events 9 and 10 in Fig. 1). During this period, earthquake occurrence in southern LV increased significantly. In contrast, the northern LV experienced a pronounced decline in seismicity, with the shallow rate dropping from 10.75 to 2.51 events/day and the deep seismicity rate also dropped from 3.79 to 2.46 events/day (Fig. S17). Stage 2 is therefore characterized by a cascade of seismic slip events occurring primarily outside the swarm zone.

Stage 3 (2023 to 2024): Gradually accelerated seismic and aseismic slip rates. Beginning in October 2022, the RESs-derived slip rate slightly increased, as shown in the long-period trends (6-month averaging windows in Fig. 6b), and continued until a few months prior to the April M7.3 event. While the annual number of repeating earthquakes remained nearly constant (~45 events/year, Fig. 6c), their seismic moment rate increased from ~1.3–1.4×10¹⁵ to 2.2 ×10¹⁵ N-m/day (Fig. 6d). This pattern indicates a larger slip per event, consistent with progressive deep fault creep and mechanical weakening at depths greater than 15 km.

Meanwhile, background seismicity from both original and declustered earthquake activities intensified, particularly at shallow depths when compared to earlier time periods. As shown by Fig. 6f–h, seismicity rate, a- and b-values increased notably in the shallow crust following the cascade of M6+ events. The rising a-value reflects an overall increase in seismic activity, while the increasing b-value suggests a greater proportion of smaller events. This pattern is consistent with Fig. S17, which shows a rise in the occurrence rate of shallow earthquakes from 2.51 to 8.37 events/day across 2023, while seismicity at greater depths shows a relatively minor increase, from 2.46 to 5.05 events/day. The declustered catalog further reveals that both shallow and deep seismicity began a gradual, slight rise roughly one year before the mainshock.

To reconcile (1) the observed increase in deep aseismic slip and RESs moment rate, (2) the rise in seismicity rate in both declustered and original earthquake catalogs, and (3) the elevated a- and b-values in the shallow crust one year before the mainshock, we propose a depth-dependent fault evolution model. At greater depths (>15 km), the CRF appears to undergo progressive aseismic slip, likely accompanied by upward fluid migration. This deep process may have weakened the overlying crust and promote more small events on both the LVF and CRF. The concurrent rise in shallow seismicity and elevated a- and b-values suggests a brittle response within a more heterogeneous and a possibly fluid-rich environment. These patterns suggest that deep slow slip and shallow brittle failure have operated as coupled processes along the CRF, which may represent an intermediate stage in the preparatory phase leading to large earthquake nucleation. However, given that aseismic slip signals are resolved only at monthly timescales, it remains difficult to establish a direct temporal link between aseismic and seismic slip in the final month preceding the mainshock.

Stage 4 (2024 April to May): M7.3 event and renewed shallow swarm activity. The April 3, 2024 M7.3 mainshock, with a hypocenter depth of 21 km, was followed by a sequence of aftershocks aligned along the steeply east-dipping LVF until April 22. At that time, a new shallow swarm was initiated at shallow depths above 11 km on the gently west-dipping CRF (Fig. S18). This post-mainshock sequence included 5 M6+ events while two of them overlapped spatially with the 2021 swarm, likely representing a recurrence of swarm activity on the same structural segment. Because the RESs catalog ends shortly after the mainshock, the evolution of deeper aseismic slip following the mainshock remains unconstrained and requires further monitoring.

In summary, we propose that since 2019, the CRF has been progressively activated through both RESs and swarm activity. As illustrated by the conceptual model in Fig. 7, the mid-2021 swarm sequence initiated at ~14 km depth, partially overlapping with the creeping segment of the CRF characterized by abundant RESs, and subsequently migrated upward to ~5 km and downward to ~ 20 km. During this period, deep aseismic slip remained elevated, suggesting that a propagating fluid pressure perturbation may have driven both the swarm and RESs activity observed in Stage 1. Significant aseismic slip likely contributed to the triggering of fluid-assisted earthquake swarms, with upward fluid migration potentially weakening the fault at shallow depths. However, the 2022 M6.1 event in Stage 2, located ~30 km south of the swarm zone and opposite the migration direction, implies that fluid migration alone cannot fully explain its occurrence. In Stage 3, seismicity intensified along a broader extent of the fault, especially at shallow depths, and was followed by the 2024 M7.3 mainshock on the east-dipping LVF and renewed swarms on the west-dipping CRF. The recurrence of CRF swarm activity after the mainshock highlights the need to evaluate how prior aseismic slip and earlier seismic events may have jointly loaded this segment. The following Coulomb stress models assess whether the observed sequence of seismic and aseismic slip episodes collectively promoted failure on the M7.3 rupture plane.

Stress triggering models
To establish whether the sequential seismic and aseismic slip episodes can be explained by static stress triggering, we applied a series of Coulomb stress models and focused on the role of aseismic slip on the M6+ events at Stages 2–4. The fault models used for individual aseismic and seismic slip events are described in Methods “Fault models for static stress computation” and listed in Table S4. The stress evolution at successive rupture sites is illustrated by sequential plots of the Coulomb stress change in Figs. S19–S22.

As shown in Fig. S19, cumulative slip from the majority of 2021 swarm sequence at depths shallower than 11 km (shallow aseismic slip episodes) resulted in a 10 kPa increase in Coulomb stress on the 2022 M6.1 rupture, while the aseismic slip from RESs at deeper portion of the fault contributed an additional 6 kPa stress change. This suggests that the M6.1 event may have been triggered by aseismic slip episodes in 2021. Both the shallow and deeper aseismic slip however, induced a negative stress change on the fault segment that later hosted the April 22, 2022 M6.7 event.

The stress changes associated with major seismic events (the 2022 March 22 M6.7, 2022 June 20 M6.1, 2022 September 18 M7.0, and 2024 April 3 M7.3 events) were also examined by evaluating the stress changes induced by preceding ruptures (Figs. S20–S22). The results indicate that the major earthquakes in 2022 can be largely explained by stress change from previous M6+ events. Specifically, the 2022 March 22 M6.7 event produced a 10 kPa stress increase on the fault segment that later hosted the 2022 June 20 M6.1 event. This M6.1 event later contributed a 2 kPa stress increase on the fault of the subsequent M7.0 event. However, the September 18 M7.0 event resulted in a −3 kPa stress change on the fault segment hosting the 2024 M7.3 event. This suggests that the M7.0 event alone does not fully account for the occurrence of the M7.3 mainshock.

Instead, the combined effect of both aseismic and seismic events likely played a crucial role in the stress triggering of the 2024 M7.3 event. This interpretation is supported by the cumulative Coulomb stress increase of up to 30 kPa at the eventual rupture location, resulting from deep aseismic episodes in 2021 and three major M6+ ruptures in 2022 (Fig. S23). Notably, a stress change on the order of 10 kPa is commonly considered sufficient to influence fault failure timing in tectonically loaded systems 40–42. Therefore, a cumulative 30 kPa increase represents a substantially elevated stress state, capable of significantly advancing the timing of a large earthquake. In addition to fault weakening due to upward fluid migration, this level of stress accumulation from both slow and fast slip along the CRF likely contributed to the enhanced seismicity and RESs activity observed during Stage 3, ultimately facilitating nucleation of the April 3, 2024 Mw 7.3 mainshock.

Methods
Earthquake catalog
We analyzed seismicity in the area bounded by latitudes 22.6°–23.4° and longitudes 120.8°–122° from January 2000 to April 2024 using data from the Central Weather Administration (CWA) in Taiwan. Crustal earthquakes with depths shallower than 50 km were selected, corresponding to the average Moho depth in eastern Taiwan (Wang et al., 2010). To ensure statistical robustness in the temporal distribution of earthquakes, we applied the maximum likelihood estimation method 72 to determine the minimum magnitude of completeness (Mc), which was estimated to be 2.0. Consequently, a minimum magnitude threshold of 2.0 was adopted for the seismicity analyses. The final catalog comprises 133,627 events with magnitudes ranging from 2.0 to 7.3.

Repeating Earthquake identification
Using the published repeating earthquake waveform procedure 30 for events from 2000 to the end of 2011, the repeating earthquake catalog was updated from 2012 to the end of April 2024. The same procedure 30,32 was followed to identify repeating earthquakes using two criteria: (1) a cross-correlation coefficient (CCC) higher than 0.9 at more than three stations and (2) a small differential S-P time (dSmP) that ensures more than 50% source overlap. By assuming 3 MPa stress drop, Vp = 4.0 km/s, and Vp/Vs = 1.78, the threshold of dSmP ≤ 0.02 s is used for M ≥ 2.5 events, while dSmP ≤ 0.01 s is used for M < 2.5 events.

The seismogram records used in this study were selected from the seismic stations operated by the Taiwan Seismological and Geophysical Data Management System (GDMS). GDMS was established through a collaboration between the Central Weather Administration (CWA) and the Institute of Earth Sciences, Academia Sinica (IESAS). Vertical component seismograms from 50,934 M≥2 earthquakes in the CWA hypocenter catalog for the period January 1, 2012 to April 30, 2024 were used. All earthquake pairs with separation distances of less than 20 km were initially selected. The time window for each event included 3 s before the P-wave arrival and 57 s afterwards.

Slip estimate for individual repeating events
Assuming that the cumulative slip of repeating earthquakes reflects tectonic loading over a single seismic cycle and can be approximated using GPS-derived interseismic slip rates 73, as previously adopted in studies 30,32,34. We combined these geodetically constrained slip rates with the cumulative seismic moment release of each sequence to estimate the average rupture area for individual repeating events, as defined in Eqs. 3 and 4. The seismic moment (Mo) was calculated using catalog magnitude (ML) data converted to moment magnitude (Mw) 74, and then to Mo using the empirical formulation 75:

MW = (0.91 ± 0.03) ML + (−0.07 ± 0.15) (1)

log Mo = 1.5 MW + 10.73 (2)

The rupture area A and seismic moment release rate (Mf0) for each repeating sequence can be computed using Eq. 3, and slip for each event in the same sequence is inferred using Eq. 4.

Mf0 = M0 / Ti = μAd / Ti (3)

d = M0 / μA (4)

Ti is the total duration of the sequence, μ is the shear modulus, and d is the interseismic slip rate. By applying this framework, a previous study established empirical scaling relationships between slip and seismic moment for different regions 30. In the southern Longitudinal Valley (LV), the repeating earthquake sequences (RESs) yielded a best-fit relation of log d = –1.21 + 0.11 log Mo, while in the northern LV, the corresponding scaling was log d = –1.96 + 0.14 log Mo. These regional differences reflect variations in fault structure, interseismic loading, and RES depth, and highlight the importance of local calibration in slip moment scaling analyses.

Earthquake swarm identification
To distinguish earthquake swarms from mainshock-aftershock sequences, the swarm sequence is expected to exhibit high seismicity density without a dominant mainshock. In more swarm-like clusters, a higher proportion of background seismicity is typically observed. A previous study 26 proposed a composite declustering method that mitigates uncertainties associated with arbitrary parameter choices in individual models. Following their approach, we applied three declustering methods: the Epidemic-Type Aftershock Sequence (ETAS) model 76, the nearest-neighbor method (NN13) 77, and the interaction-zone method (Re85) 78, which links earthquakes based on spatial and temporal proximity.

The identification of earthquake swarms requires the determination of density rate anomalies (δμ0). To do so, the density rate (μ0) of the declustered seismicity is first estimated to represent the background rate as

μ0(x,y,t) = Σ ωi e−((x−xi)2 + (y−yi)2) / 2πl2 e−|t−Ti| / T × 1 / 2αi (5)

where ω denotes the earthquake probability of the background event, l denotes a spatial smoothing parameter, and T indicates the temporal smoothing parameter. ti represents the occurrence time of an earthquake in the catalog. The scaling parameter αi is defined as

αi = 2T (1 − e−(ts−ti)/T) (1 − e−(te−ti)/T) (6)

where ts and te are the starting and ending times of the input catalog. The procedure was followed for estimating the stacked density rate (μ0stack) using three different declustering algorithms 26. The same parameter settings for the declustering algorithms were maintained, but the smoothing parameters for the density rate (Eq. 1) were modified from T = 15 days and l = 30 km to T = 5 days and l = 3 km, in consideration of the different magnitude cutoffs. The chosen cutoff magnitude was 3.0 for all of Taiwan 26, but here M 2.0 is used to obtain a more complete swarm catalog in eastern Taiwan.

To identify density rate anomalies (δμ0), the local rate change of μ0 in each swarm was estimated as δ(t) = {(μ0stack(t) - μ0average(t)) / μ0average(t)} x 100, where μ0average(t) is the averaged μ0(x,y,t) within a 3 km radius of each event, with the average μ0 being computed over a two-year window centered on time t (Fig. S26a). As shown by Fig. S26b, the resulting number of clusters is found to decrease smoothly with increasing threshold on δ. To discriminate between swarm-like and aftershock-like clusters, the latter type is expected to have greater magnitudes. Thus, in Fig. S26c, the number of clusters with a maximum magnitude greater than 6 is plotted against δ threshold. The number of M≥6 clusters decreases from 13 to 0 when δ changes from −20 to 15, while the corresponding number of total clusters decreases from 140 to 27. To retain more swarm-like clusters while excluding all the aftershock-like clusters, the optimal δ was chosen as 15 in this study.

GNSS data processing and analysis procedures
The Taiwan Continuous GNSS Array is operated by the Institute of Earth Sciences, Academia Sinica (IESAS), the Central Weather Administration (CWA), the Geological Survey and Mining Management Agency (GSMMA), and the Ministry of the Interior (MOI). GNSS position time series for stations across Taiwan are available for download from the Taiwan Geodetic Model (TGM) website (https://tgm.earth.sinica.edu.tw/TseriesList). GNSS data are processed using GipsyX/RTGx software 79, incorporating final orbit and clock parameters from Caltech Jet Propulsion Laboratory (JPL) for precise satellite positioning. To mitigate atmospheric errors, we apply an ionosphere-free linear combination technique 80, which integrates multiple frequency observations to reduce ionospheric effects. Tropospheric delays are corrected using the Vienna Mapping Function 81. Ocean tide loading effects are removed using the FES2004 model 82, which accounts for deformations caused by Earth-ocean gravitational interactions. The daily coordinate time series are processed in the International Terrestrial Reference Frame 2014 (ITRF14) and then transformed into a local NEU (north, east, up) coordinate system.

To analyze these GNSS timeseries, we first apply least-squares regression to model and remove long-term linear trends, annual and semi-annual seasonal signals, instrumental offsets, and coseismic displacements. For earthquakes showing clear co-seismic offsets, we further correct for postseismic deformation by fitting an exponential decay function of the form:

Δ(t) = Ae−t/τ (7)

where Δ(t) presents the time-dependent postseismic displacement, A is the initial postseismic amplitude, and τ is the characteristic decay time. Only M≥6.0 events in this study area (listed in Tables S1 and S5) are assumed to generate postseismic deformation considered in this correction. For smaller events, we evaluate the significance of potential postseismic effects by comparing the root-mean-square (RMS) misfits between models with and without exponential decay term. If the difference is less than 10% of the mean observational uncertainty, we consider the postseismic effect negligible and did not apply postseismic correction. The fitting of the exponential decay capture the postseismic signal observed in the time series at each station. The estimated postseismic decay constants (τ) and their uncertainties are summarized in Table S5. Residual time series are computed by subtracting the full regression model (including postseismic terms) from the raw observations. The original and corrected GNSS time series are shown in Fig. S14 and Fig. 5, respectively.

Fault models for static stress computation
To better understand whether the sequential seismic and aseismic slip episodes can also be explained by static stress triggering, a series of computations of Coulomb stress change was conducted. The source model and receiver fault assumptions are listed in Table S4. The apparent 2021 slow slip episodes comprised two segments: one illuminated by the shallow swarm activity (5–15 km) and the other by the deeper RESs (15–25 km). The fault geometry along the northern CRF was approximated by the averaged strike, dip, and rake from the M5+ events, while the fault length and width were approximated by the spatial distribution of both the swarm and RES activity. Similarly, the aftershock distribution and published source model were adopted for parameterizing the three mainshocks in 2022. The fault areas of the M6.1 June 20, 2022, M6.7 March 22, 2022, and M7.0 September 18, 2022 events are represented by the spatial distribution of the aftershocks that occurred in the month following each mainshock, while their fault orientations were derived from the focal mechanisms of the mainshocks.

Data availability
All data are available in the main text or the supplementary materials. The Python software package Obspy (www.obspy.org) was used for seismic data processing and waveform filtering. The stastic Coulomb stress change is calculated using Coulomb 3.4. Figures were produced using Generic Mapping Tools (GMT), and Matlab. The catalogs of repeating earthquake sequences, earthquake swarms are available at https://doi.org/10.5281/zenodo.14064472. The earthquake catalog of Taiwan can be retrieved from https://gdmsn.cwb.gov.tw/. The focal mechanism of Taiwan can be retrieved from https://bats.earth.sinica.edu.tw/.

Code availability
All software used in this work is open source. The code used to generate each figure and result is available through the contact information from the original publications. Requests for further materials should be directed to Wei Peng (weipeng@ntnu.edu.tw).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Peng et al. study the evolution of microseismicity leading up to the M7.3 Hualein city earthquake with a focus on repeating earthquakes and swarms. They use a 24 year long catalog along a large portion of the longitudinal valley fault. The analysis of the catalog centers on inferring the presence of aseismic slip from repeating earthquake sequences and seismic swarms. The abstract claims an increase in aseismic slip rate for ~3 years prior to the M7.3 earthquake which would be a very important result. The primary argument for this would appear to come from a combination of figure 3b (repeating earthquakes) and figure 4 b (swarms). It is very hard in both figures to see the exact behavior leading up to the M7.3 earthquake. Similarly, the black curve (RES) in figure 6a does not show any anomaly before the M7.3 earthquake? My primary reaction to this paper is that after lots of zooming in and out of the figures it is not at all obvious to me what in the analysis supports what would be such a widely important conclusion. Perhaps it’s a matter of presentation, as the 3-year anomaly is not really highlighted on any figures, etc., but it’s really not a convincing paper in its current form. It seems they have clear evidence of afterslip following some of the M6 earthquakes but it’s not at all clear where the claim of an ‘acceleration of aseismic slip rate’ comes from. The lack of clarity in the figures combined with the complete lack of validation of the aseismic slip (for or against) with geodetic data make this a very incomplete paper. I recommend that it would require major revisions with a focus on clearly documenting the primary conclusion.\nFigure 4b is presented as the key argument that increased aseismic slip is involved in/preceeds most of the earthquake swarms in the catalog (line ~200). Swarms like these no doubt identify time windows when something interesting is happening at depth and fluid and aseismic slip transients are the most likely explanations. I’m not really sure what the authors see in Figure 4b that leads them to their conclusions. It seems that only two of the swarms would be candidates for having aseismic slip during the swarm and most show no anomaly either before or after the swarm. Mostly Fibure 4b would seem to show a lot of variation. Similarly, Figure 4c would seem to indicate that most of the large quakes are not associated with any swarm activity preceeding them. 4c could use a zoomed in view of 2024 similar to figure 2b. Hard to tell what happens relative to the M7.3 mainshock. Similarly Figure 5c doesn’t seem to show any clear expansion over time or fit the diffusion curves. Overall, it’s not clear that the swarms in Figures 4 and 5 are indicative of a clear increase in aseismic slip.\nI’ve never fully bought into the simplistic mapping between repeating earthquakes and aseismic creep rates over a wide portion of the same fault for a number of reaons not the least of which is the very localized nature of M2s. These types of highly localized observations are often mapped into generalized statements about entire fault segments. For instance, Materna, Taira, and Burgmann [2018] used the same methodology as this study to extrapolate from a narrow strip of repeating earthquakes to conclude that large portions of the Mendocino transform fault were failing primarily by aseismic creep. Yet this portion of that fault just ruptured in a M7 earthquake a few years later clearly having built up slip deficit over the preceding decades when it was claimed to be creeping.\nIt is very strange that this is a paper that is fundamentally trying to infer aseismic slip but does not evaluate any geodetic data or even compare to detailed studies of the same fault. For instance, the Thomas et al. 2014 JGR paper has extremely detailed maps of locking and aseismic slip along these exact fault segments covering over half the study period. Why are there no comparisons (or even mention) to those aseismic slip rates? I’m not an expert in the LVF, but I’m sure there are more recent analyses of the very dense geodetic networks there that could inform the later part of the study period.\nThe section from lines 81-130 is hard to follow and quite detailed. I’m not sure if there’s a better way to do this, but it amounts to a distraction from the key points of the paper.\nThe key relations on line 147 make very specific inferences about the stress drop distributions on the two faults. These should be spelled out explicitly and justified as they are the underpinning of the primary calculation in the paper. This is also buried in the ‘dSmP’ criteria for identifying repeaters. Whether they have 50% source overlap is very much a stress drop question. All of this should be more explicitly stated in terms of the underlying assumptions about stress drop with magnitude, etc.", "reviewer_2": "The manuscript written by Peng, Chen, and Burgmann describes their observations of the interaction of repeating earthquake-derived aseismic slip and seismic swarms leading to a large M 7.3 earthquake in eastern Taiwan. The overall presentation is clear and easy to understand. The data quality and presentations are highly reliable, providing a robust observation for the precursory slow slip before a large earthquake. This manuscript is significant for studying precursory slip leading to a large earthquake.\nBased on the observations, this manuscript discussed the possible controlling factors for earthquake swarms' diffusivity and fluid's origin. However, I do have one suggestion for the manuscript. Considering the 24-year catalog of repeating earthquakes, swarms, and regular earthquakes, at least five episodes of moderate to large earthquakes occurred in eastern Taiwan. With this dataset, a few parameters can be derived relating to their locations and timing: the aseismic slip at depth calculated from the repeating earthquake sequences, earthquake swarm diffusivity, and the occurrences of moderate to large earthquakes. The authors can derive an empirical relationship among these parameters and see how they relate. This relationship or formula may explain how the precursory slip leads to a large earthquake.", "reviewer_3": "This paper attempts to gain insight into aseismic processes related to the occurrence of Magnitude 6-7 class earthquakes by extracting repeating earthquakes and earthquake swarm activities from long-term seismic observation data.\nThe authors claim to 'observe an acceleration of aseismic slip rate that occurred approximately three years preceding the April 3, 2024 Mw 7.3 Hualien earthquake'. While the proposed correlation between major earthquakes and aseismic slip is potentially very interesting, the supporting evidence presented in the current form falls short of being conclusive.\nThe main text is written based on the assumption that both repeating earthquake sequences (RES) and earthquake swarms are driven by aseismic slip (or slow slip), manifesting as seismic activity as RES and earthquake swarms. However, in general, the relationship between earthquake swarms and aseismic slip is not obvious and should be supported by additional observational data, such as geodetic measurements. Such data is not presented in the paper. Moreover, given the mention of the involvement of fluids in earthquake swarms, it is conceivable that fluid migration alone could explain earthquake swarm activity. The necessity of invoking aseismic slip in addition to fluid migration has not been adequately justified.\nThe discussion of hypocenter migration during swarm episodes is dependent on the definition and extraction methods of earthquake swarms. While Figure 4c plots only the extracted swarm activities, it is difficult for readers to judge the validity of the extraction within the overall seismic activity. For example, it would be necessary to show a magnitude vs. time plot similar to Fig. S2b within the rectangular area outlined in gray in Fig. 4a.\nEven if we accept the extracted earthquake swarms, it is desirable to present spatiotemporal plots where the hypocenter locations are projected onto axes representing north, east, depth, or the overall migration direction, in addition to plotting the distances between earthquakes as in Fig. 5c. In the case of Fig. 5c, the main text states that the hypocenters migrate from the deeper to the shallower part of the CRF. However, is it appropriate to discuss the migration using plots with axes that do not necessarily reflect such a directional movement?\nFurthermore, the necessity of applying the diffusion model is not clear. For example, in Fig. 5c, the location of the initial event is unclear (it should be clearly indicated in the figure which earthquake is the reference), and the events have already spread over 20-30 km within the first 10-20 days, but this part seems to be ignored. Additionally, in Fig. 5b, the initial earthquakes appear to occur in the middle (10-17 km) of the overall hypocenter distribution. If so, it suggests that the migration is not only towards a shallower part but also towards a deeper part, but the latter is not discussed.\nFurthermore, even if the diffusion model can be used for explaining the migration of the hypocenters in the extracted earthquake swarm activities, the diffusion model fitting shown in Figures 5c and S5-S13 is performed with respect to the distance from the first EQ of an extracted earthquake swarm cluster. Therefore, if a different earthquake within the same earthquake swarm activity had been chosen as the initial EQ, the spatial location of the reference for \"distance\" could change significantly, potentially leading to a large difference in the estimated value of D. From this perspective, as mentioned above, it is necessary to consider the examination along a distance axis in one direction within the 3D space.\nRegarding the case of Figure 5c, it could be argued that the sequence of seismic activity can be divided into three approximately distinct phases: an initial period of approximately 20 days, followed by a period of about 10 days starting around the 50th day, and a final period of about 40 days commencing around the 70th day. Even within each of these phases, hypocenter migration can be observed. It is therefore worth considering whether the possibility of such migration within these sub-clusters should be further investigated.\nminor points:\nThe authors set a minimum magnitude threshold of 2.0 for the earthquake catalog. A more in-depth discussion is warranted regarding the implications of this setting (Mc) on the analysis of earthquake occurrence rates. It is recommended that this topic be explicitly addressed in the Methods section.\nFigure 1:\nThe focal mechanisms (beach balls), along with the text of earthquake year and Mw, are overlapping with symbols such as open circles and brown diamonds, making these symbols very difficult to see. The focal mechanisms and related text should be moved to a position further out, where they do not overlap with other symbols.\nIn Figs. 3 and 4, the \"along strike\" distance is shown, but it is unclear where this corresponds to within Figure 1a. A distance axis should be added to Figure 1a to clarify this.\nSince there are several different symbols used in the figure, a legend should be included to explain their meanings.\nThe black rectangular regions and the square bracket-like shape to the right of the letters A' and B' are not explained in the figure caption.\nSimilarly, in Figure 1b and 1c, the symbols are difficult to see due to overlapping text. There is also overlapping text related to different earthquakes, which is meaningless.\nGiven that Table S1 lists major earthquakes, it would be helpful to add numbers to these events and only show these numbers in the figure, or to use some other method to distinguish between different earthquakes.\nL82-89, Figure 2:\nIt is unclear which range of earthquakes within the CRF and LVF directions (generally north-south) are sampled for the plot. Please indicate in Figure 1 the specific area for which the earthquake density is calculated.\nThe density is described as \"each 1 square kilometer area and within 20-day time,\" but the color scale only uses \"low\" and \"high\" without specific numerical values. Please provide specific numerical values for the color scale.\nThe labels for major earthquakes (year, Mw) overlap with the density colors, especially for earthquakes after 2022, making it difficult to see the density. Please improve the visibility by repositioning the labels.\nIn L89, \"vertical red bar\" is mentioned, but it is unclear which part is being referred to. Additionally, there are several vertical blue dotted lines in the figure, but there is no explanation for these in the caption.\nTo demonstrate the simultaneous activation of seismic activity on both sides of the fault, I would suggest that the authors add a plot that shows the time-dependent cumulative number of earthquakes within defined areas on either side of the fault in a map similar to Figure 1.\nIn addition to the overall 25-year period, zoomed-in plots for each of the five specific periods of interest should be included in the supplementary information for both Figure 2a (density plot) and the suggested cumulative number of earthquakes plot.\nL85:\nThe abbreviation \"LV\" is used before its definition is given in L115.\nL101-104, Figure S1:\nAlthough the authors discuss aftershock activity in the LVF and CRF regions, while the cumulative seismic moment they present is interesting, it would be more common and informative to also show the time series of the cumulative number of earthquakes and magnitude-time plots. Additionally, in Figure S1, only the long-term changes over 25 years from 2000 are plotted, but for a more detailed evaluation of aftershock activity, it would be more appropriate to show zoomed-in plots for shorter time periods.\nAs pointed out for other figures, the presence of symbols like A-F and labels indicating the occurrence year of major earthquakes can obscure the plotted quantity (cumulative seismic moment). In particular, the \"D\" symbol obscures a line, which is unacceptable.\nMovie 1:\nI initially assumed that the color of the symbols represented the occurrence time of earthquakes. However, the presence of additional colors like gray and black makes it difficult to interpret the data. Furthermore, the limited colors for the circle symbols, primarily pink and blue, making it difficult to distinguish the differences in occurrence times. The overall purpose of this movie is unclear.\nMoreover, despite being a movie format, the lack of time or date labels for each frame makes it very difficult to understand what the author wants to show. I think it would be more clear if the time scale in the lower left corner had an arrow indicating the corresponding time, and the arrow moved to the right as time progresses.\nFigure 3a:\nAs I comment on Figure 1, the text labels indicating year and Mw of earthquakes overlaps with other symbols, making them very difficult to see.\nThe figure caption mentions \"shaded and open polygons,\" however, it is not clear which specific features of the figure these terms refer to.\nFigure 3b:\nI believe that overlaying a color map of earthquake density with the RES slip rate would compromise the clarity of the figure. I recommend presenting these two datasets separately.\nFurthermore, the vertical axis of the RES slip rate plot should be clearly labeled, particularly for the creeping CRF slip rate plotted at the top. The current plot makes it difficult to identify the zero level (0 cm/day), and the data appears to extend into negative values. Additionally, the labels indicating major earthquakes overlap with the data lines, especially after 2022, making it difficult to see. These issues should be addressed before finalizing the figure.\nThe color scale in this figure is described as representing \"event density\" in the figure caption. However, given the presence of the RES slip rate within the same figure, this could be interpreted as referring to the event density of repeating earthquakes. To avoid any ambiguity, please provide a more specific description that clearly distinguishes between these two quantities.\nFigure S2b:\nThe cumulative number of EQs and the magnitude-time plot overlap, and the labels indicating the years and Mw of major earthquakes, as well as the colored sections for different time periods, further obscure the background magnitude-time plot. As mentioned in other figures, please improve visibility by plotting these elements separately.\nL156-159:\nWhile the author states a correlation between the aseismic slip rate and RES occurrence rate around 2003-2004, a direct comparison of Figure 3b and Figure S3b is required to verify this claim. However, both figures are plotted over very long periods, making it difficult to compare the changes in the 2-3 year period mentioned here. I would recommend creating a new figure with Figure 3b and Figure S3b plotted vertically side-by-side, sharing a common x-axis (time axis) and zooming in on the 2-3 year period of interest.\nFigure S3:\nThe M-T plot and the cumulative RE number line overlap, making it very difficult to see. To improve clarity, I suggest dividing the figure into two panels, sharing a common x-axis.\nWhile the figure caption mentions \"The green, red, blue sections,\" the use of color alone to represent different time periods may be difficult to distinguish. I think it would be better to use a shape like a square bracket (]) at the top of the top panel to indicate the time periods.\nAdditionally, please provide a detailed explanation in the figure caption for all elements, including arrows, red dashed lines, and numerical values.\nL199-200:\nThe authors excluded swarm sequence #10 from their discussion because they have not updated the RES catalog.\nHowever, I think it would be better to update the catalog and include the swarm sequence #10 in the analysis. The current approach seems rather incomplete.\nL200-201:\nThe authors states that \"6 out of 9 swarm sequences were preceded by the increasing aseismic slip rate,\" but Figure 4b shows that only #1, #2, #3, and #8 exhibit this pattern. For sequences #7 and #9, the aseismic slip rate and the swarm onset appear to occur simultaneously.\nAdditionally, the confidence interval of the RES slip rate is also unclear.\nFurthermore, the selection of swarms is dependent on the definition used in this study. The robustness of such arguments should be examined if the definition is changed.\nFigure 4:\n(a) Similar to the other figures, the text overlaps with the symbols.\n(c) The overlay of swarm locations (along strike distance) and the RES slip rate results in portions of the black line being obscured. To improve clarity, the plotted swarm locations should be adjusted vertically to avoid overlapping with the black line.\nIt should be clearly stated in the caption whether the plotted RES slip rate (black line) is identical to one of the two RES slip rates plotted in Figure 3b, or if it represents a different value. The current state of the figure is likely to confuse readers.\nThe diameters of the open circles vary, but there is no explanation for what these differences represent (similarly in Figure 1).\nMovie 2:\nAs pointed out in Movie 1, it would be better to show the time of each frame in the movie.\nSimilarly, the color scheme for the symbols should be improved as I pointed out in Movie 1.\nThe sudden appearance of numerous points at the beginning of the movie creates a strange impression as it is unclear how simultaneous they are because the time interval between each frame is unclear.\nL219:\nThe text mentions \"low Vp and Vs,\" but this data is not presented in any of the figures. At least a citation should be provided.\nL277: \"located to the south of the earthquake swarms, opposite to their propagation direction\"\nPrior to this sentence, the text has described the vertical migration of hypocenters from depth to a shallower part. The abrupt introduction of the term \"south\" seems to imply a migration in north-south direction. This ambiguity makes it difficult to determine the exact direction of migration.\nFigure 6a:\nAlthough four quantities are plotted in this figure, it is very difficult to understand due to the unclear vertical axes. To improve clarity, I suggest dividing the data into four separate subplots, each with a clear vertical axis and a shared horizontal time axis.\nFigure S4:\nSame comments as for Figure 1.\nFigures S5-S13:\nSame comments as for Figure 5.", "reviewer_author_response": "We sincerely thank the reviewers and handling editor for their thoughtful and constructive comments, which have significantly improved the quality and clarity of our manuscript. In this revised version, we have addressed all reviewer concerns with detailed point-by-point responses and have made substantial improvements to the text, figures, and supporting analyses. Key revisions include expanded GNSS-based deformation analysis, more rigorous statistical testing of swarm and aseismic slip interactions, improved figure clarity and consistency, and a clearer explanation of the methodologies and interpretations throughout the manuscript. Below, we provide a detailed response to each comment, with line numbers corresponding to the revised manuscript. Changes made in response to reviewer suggestions are highlighted in grey (track change file) to facilitate evaluation.\nResponses to the comments by Reviewer 1\nQ 1-1 Peng et al. study the evolution of microseismicity leading up to the M7.3 Hualein city earthquake with a focus on repeating earthquakes and swarms. They use a 24 year long catalog along a large portion of the longitudinal valley fault. The analysis of the catalog centers on inferring the presence of aseismic slip from repeating earthquake sequences and seismic swarms. The abstract claims an increase in aseismic slip rate for ~3 years prior to the M7.3 earthquake which would be a very important result. The primary argument for this would appear to come from a combination of figure 3b (repeating earthquakes) and figure 4 b (swarms). It is very hard in both figures to see the exact behavior leading up to the M7.3 earthquake. Similarly, the black curve (RES) in figure 6a does not show any anomaly before the M7.3 earthquake? My primary reaction to this paper is that after lots of zooming in and out of the figures it is not at all obvious to me what in the analysis supports what would be such a widely important conclusion. Perhaps it’s a matter of presentation, as the 3-year anomaly is not really highlighted on any figures, etc., but it’s really not a convincing paper in its current form. It seems they have clear evidence of afterslip following some of the M6 earthquakes but it’s not at all clear where the claim of an ‘acceleration of aseismic slip rate’ comes from. The lack of clarity in the figures combined with the complete lack of validation of the aseismic slip (for or against) with geodetic data make this a very incomplete paper. I recommend that it would require major revisions with a focus on clearly documenting the primary conclusion.\nAseismic slip and seismic swarms leading up to the 2024 M7.3 Hualien earthquake\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nResponse 1-1\nThe reviewer raised concerns regarding (1) the lack of validation of the acceleration of aseismic slip rate using geodetic measurements and (2) the level of clarity in the presented evidence and figures documenting the seismicity rate changes and the potential role of postseismic afterslip in producing the observed features. For this revision, we have incorporated GNSS data analysis (new Section “4. Geodetic observations during the swarm sequences”) to evaluate for aseismic slip signals. Please see the processed GNSS time series in the new Figure 5 and corresponding lines 212-256 for details.\nWe also revised the depiction of the aseismic slip behavior in the 24-yr period, and we extended the highlighted period of our analysis from “2021-2024” (based on the onset of the longest swarm sequence) to “2019-2024” (based on the significant increase of aseismic rate during this time interval). As shown in Figure R1 below and new Figure 2a, the RES-derived aseismic slip rate is found to have accelerated starting from 2019, coinciding with the time of 2019 Mw6.1 event. The key features of seismic and aseismic activities summarized in Figure 6 are now shown over 5 years. Please see the text in lines 259-333 for details.\nFigure R1 (new Figure 2a). Magnitude distribution (top) and occurrence rate (bottom)\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nof RESs in the area of northern Longitudinal Valley (NLV). Vertical dashed lines indicate the occurrence time of Mw ≥ 6 events. Black arrows in the bottom panel denote the significant change of the RESs occurrence rate, breaking the 24-yr into three periods with different averaged occurrence rate as denoted by light green, blue, and red areas, corresponding to 41.7, 33.3, and 45.8 events/yr, respectively.\nQ 1-2 Figure 4b is presented as the key argument that increased aseismic slip is involved in/preceeds most of the earthquake swarms in the catalog (line ~200). Swarms like these no doubt identify time windows when something interesting is happening at depth and fluid and aseismic slip transients are the most likely explanations. I’m not really sure what the authors see in Figure 4b that leads them to their conclusions. It seems that only two of the swarms would be candidates for having aseismic slip during the swarm and most show no anomaly either before or after the swarm. Mostly Figure 4b would seem to show a lot of variation. Similarly, Figure 4c would seem to indicate that most of the large quakes are not associated with any swarm activity preceeding them. 4c could use a zoomed in view of 2024 similar to figure 2b. Hard to tell what happens relative to the M7.3 mainshock. Similarly Figure 5c doesn’t seem to show any clear expansion over time or fit the diffusion curves. Overall, it’s not clear that the swarms in Figures 4 and 5 are indicative of a clear increase in aseismic slip.\nResponse 1-2\nWe appreciate the reviewer’s concerns regarding the acceleration of aseismic slip rate before/after the swarm sequence. In the revised manuscript, we have provided a new sensitivity test for a better presentation, as detailed below.\nIn the old Figure 4b, we tried to understand if there exists an acceleration in slip rate (from RESs) before the onset of the swarm sequences, and how close in time such acceleration occurred from the swarm. We agree with the reviewer that the previous presentation was somewhat subjective. Now we have conducted a sensitivity test, where the change of RES slip rate is obtained between a time window before (d ) pre and after (d ) the swarm onset, calculating the change of slip rate as df = d - d . post post pre\nThis df is plotted using various time windows, ranging from 1 to 12 months in Figure R2 (same as the new Figure S12). Positive values indicate that a slip rate acceleration occurred after the onset of the swarm sequence, while negative df values indicate a slip rate acceleration occurred prior to the swarm sequence. This quantitative test suggests that the description in the original manuscript of “6 out of 9 swarm sequences were preceded by increasing aseismic slip rate”, in fact, remains true for time windows of 3 – 6 months from the onset of the swarm sequence. Please see lines 186-197 in the revised text.\nFigure R2 (new Figure S12). Sensitivity test of slip rate change measured over different time windows. Here, slip rate is averaged in time windows before (d ) and after (d ) pre post\nthe swarm onset. Different colors represent different repeater sequences, with the seq. ID indicated by the number on the top. Number labels on the left indicate the number of sequences characterized by higher slip rates before the swarm onset.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nQ 1-3 I’ve never fully bought into the simplistic mapping between repeating earthquakes and aseismic creep rates over a wide portion of the same fault for a number of reaons not the least of which is the very localized nature of M2s. These types of highly localized observations are often mapped into generalized statements about entire fault segments. For instance, Materna, Taira, and Burgmann [2018] used the same methodology as this study to extrapolate from a narrow strip of repeating earthquakes to conclude that large portions of the Mendocino transform fault were failing primarily by aseismic creep. Yet this portion of that fault just ruptured in a M7 earthquake a few years later clearly having built up slip deficit over the preceding decades when it was claimed to be creeping.\nResponse 1-3\nThe conceptual model in old Figure 6b (now Figure 6f) meant to capture the along-dip distribution of swarm sequences and RESs, to compare with the observations in Figure 6a to Figure6e. It does not attempt to provide a detailed view of aseismic vs. seismic slip on the fault, which requires more comprehensive understanding of the interseismic slip model. In reality, the fault geometry is more complex in order to host both swarms and large earthquakes on the CRF. Particularly, at shallow depths particularly, the fault structures likely evolve progressively, accommodating both large earthquakes and active swarm sequences. This is now clarified in the text associated with the interpretation of Figure 6f (line 282-333).\nQ 1-4 It is very strange that this is a paper that is fundamentally trying to infer aseismic slip but does not evaluate any geodetic data or even compare to detailed studies of the same fault. For instance, the Thomas et al. 2014 JGR paper has extremely detailed maps of locking and aseismic slip along these exact fault segments covering over half the study period. Why are there no comparisons (or even mention) to those aseismic slip rates? I’m not an expert in the LVF, but I’m sure there are more recent analyses of the very dense geodetic networks there that could inform the later part of the study period.\nResponse 1-4\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\n1. In Thomas et al. (2014), the interseismic slip model was inferred for the east-dipping LVF. This study, however, is focused on the west-dipping CRF. Unlike the LVF, the CRF is poorly constrained due to the lack of geodetic data coverage in the mountain area (the hanging wall of the CRF). As shown in Figure R3, the inferred inter-seismic slip rate in the northern LV area is relatively low (Figures R 3b-c) while the RESs are located at depths greater than 10 km (Figure R3a). Our earlier work demonstrated that the contribution of the CRF affects the kinematic modeling of the LVF (Peng et al., 2023), and the uncertainties in fault parameters chosen for the CRF are large, requiring better constraints from relocated seismicity and better coverage of geodetic measurements in the eastern flank of the Central Range. What we have added to the current study, is an exploration of the aseismic slip episodes using the existing GNSS data. This is demonstrated in the new section on “4. Geodetic observation during the swarm sequences” (see lines 212-256 and Figure 5) and detailed below (after Figure R3).\n[REDACTED]\nFigure R3 (not included in the manuscript). (a) Distribution of RESs (brown diamonds) in eastern Taiwan. The RESs in the southern LV area occurred along the east-dipping LVF, while in the northern LV area they mainly occurred along the west-dipping CRF.\n(b-c) GPS and InSAR data used for the inter-seismic model in Thomas et al. (2014). Note that in the northern LV the velocity gradient across the LV are low compared with the velocities in the south.\n2. The 2021 swarm sequence includes 9 M≥5 events, the concentration of swarm activity is denoted by the vertical lines in Figure R4, coinciding with the time of two largest events (M5.9 and M5.5). The raw GNSS data in Figure R4 below shows\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nthat at the time of the M5.9 even, co-seismic offsets are significant at nearby stations (WARO, FLNM, and DNFU). As highlighted by black arrows, 15 mm southward and 10 mm eastward displacement is observed at WARO station. At FLNM and DNFU southward displacements of 20 mm and 10 mm are observed, respectively. The southeastward displacements are consistent with the left-lateral focal mechanism with a moderate thrust component. As the surge of swarm activity moved to the northeast towards shallow depth, a sudden change in displacement is observed at the time of M5.5 event; Clear displacements can be found at stations TUNM, YENL, NDHU, and SOFN, with southward and eastward displacements of ~10 mm. During the earthquake swarm, additional gradual motions are obvious in the east component at YENL, NDHU, and SOFN stations, which may be related to aseismic slip and/or afterslip of the M5.9 event.\nWhen the long-term trends, annual cycles, and co-seismic offsets are removed, the corrected data (Figure R14 and new Figure 5) show visible changes spanning the 2019 swarm sequence at stations FLNM and DNFU and at stations TUNM, NDHU, and YENL during the 2021 swarm sequence. The southeastward displacements are observed at the stations close to the swarm distribution, suggesting that they are of tectonic origin instead of noise.\nFigure R4 (Figure S14). (a) Distribution of the selected 11 GNSS stations (triangles) and the swarm sequences in 2019 (red dots) and 2021 (black dots). Stars indicate the Mw ≥ 5.5 events since 2019. (b-c) (Top) N and E components of the GNSS original signals, color coded by different stations shown in (a). Vertical lines denote the\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\noccurrence time of the Mw ≥ 5.5 event since 2019. (Bottom) Magnitude-time plot with cumulative number (Purple line) of swarm events. The red, black, and grey dots in the bottom panel indicate the 2019 swarm sequence, 2021 swarm sequence, and background seismicity, respectively.\nFigure R5 (new Figure 5). Same as Fig. R4 but the long term trend, annual cycle, and co-seismic offsets were removed from the original GNSS data.\nQ 1-5 The section from lines 81-130 is hard to follow and quite detailed. I’m not sure if there’s a better way to do this, but it amounts to a distraction from the key points of the paper.\nResponse 1-5\nThis section described the seismicity on both the CRF and LVF, providing basic understanding of the seismic behavior in the study area. However, we agree with the reviewer that the details addressed here may not be relevant to the following sections. We have decided to only keep the general introduction of major earthquakes and removed the content related to details on LVF-CRF interactions from the main text.\nQ 1-6 The key relations on line 147 make very specific inferences about the stress drop distributions on the two faults. These should be spelled out explicitly and justified as they are the underpinning of the primary calculation in the paper. This is also buried\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nin the ‘dSmP’ criteria for identifying repeaters. Whether they have 50% source overlap is very much a stress drop question. All of this should be more explicitly stated in terms of the underlying assumptions about stress drop with magnitude, etc.\nResponse 1-6\nThe reviewer raised a question regarding the role of stress drop in the slip determination (using logd = − 1.21+0.11 logMo and logd = − 1.96+0.14 logMo) and identification of repeaters using the dSmP threshold. These issues were discussed at some length in Chen et al. (2020). Below we provide some details to answer this question.\n(1) What is the role of stress drop in the RESs slip determination?\nThe slip estimate for individual repeaters was inferred without making assumptions about the stress drop. Following Nadeau and Johnson (1998), we assumed the cumulative slip is a result of tectonic loading from the surroundings over one seismic cycle and can be represented by the GPS-derived interseismic slip rate in our earlier study in Chen et al. (2020). Therefore the GPS-derived interseismic slip rate together with the cumulative moment release smoothed over the duration of one seismic cycle are used to approximate the rupture dimension of individual repeater. The rupture radii derived from the seismic moment rate and the assumed geodetical long-term slip rate in both equations of logd = − 1.21+0.11 logMo and logd = − 1.96+0.14 logMo are referred to as 𝑅 (how we obtained these two equations is now introduced in Method, !\"̇\n-Slip estimate for individual repeating event). The question now is, how is 𝑅\n!\"̇\ncompared with the rupture radius assuming circular rupture from Kanamori and Anderson (1975)?\nAssuming 3 MPa stress drop and used the equation below, we are able to obtain the slip as well as 𝑅 .\n∆%\n7𝜋𝑀\n&\n𝑅 =\n∆% 16∆𝜎\nFrom the calculation in Chen et al. (2020), 𝑅 is found to be larger than 𝑅 , as their\n∆% !\"̇\nratio (𝑅 /𝑅 ) ranges from 1 to 10. That means in order to produce the rupture\n∆% !\"̇\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\ndimension considered in the slip vs. Mo scaling for RESs, 10~30 Mpa stress drop may be needed for all events. This is consistent with the averaged stress drop of 6.7-19.1 Mpa determined from empirical Green’s function using two M4 repeating sequences in eastern Taiwan (Chen et al., 2016).\n(2) How the stress drop influences the dSmP threshold for repeater identification?\nThe dSmP threshold is sensitive to (1) the magnitude of target events (2) assumption of stress drop, as addressed in Chen et al. (2020). The source size determined using\n∆𝜎 =\n'(!\nfor different magnitude gives 57-566 m for M2-4 events assuming 3 Mpa\n)*+\"\nstress drop. A stress drop of 100 Mpa, however, would lead to rupture radii of 18-176 m. To reach for more precise stress drop assumption, we consider the stress drop estimates based on empirical Green’s function for two M4 repeating sequences in eastern Taiwan (Chen et al., 2016). In their work, a peak stress drop of 64.4-155.1 Mpa was obtained for a quasi-periodic sequence with an averaged stress drop of 6.7-19.1 MPa. When the dL (separation distance between the centroids of two similar-sized events) is equal to the rupture radius r, the percentage of source overlap is only 40%, assuming circular rupture for two equally sized events. To ensure at least 50% source overlap for the repeating events, the normalized distance (i.e., the inter-event distance dL divided by the sum of the rupture radii of two events) should be smaller than 0.4, as denoted by horizontal dashed lines in Figure R6. By asking for 50% source overlap, stress drops of 3 MPa and 19 MPa lead to the thresholds of dSmP ≤ 0.012 s and dSmP ≤ 0.006 s for M2.0 events. However, the sampling rate of seismograms is 100 Hz, leading to lower resolution in time if we push below 0.01 s. Therefore, we adopted the relatively loose threshold on the rupture dimension (i.e., a 50% source overlap criteria), using 3 MPa stress drop. Consequently, in this study we adopted the method proposed in Chen et al. (2020) by assuming 3 MPa stress drop, Vp = 4.0 km/s, and Vp/Vs = 1.78. For M ≥ 2.5 and M < 2.5 events, the thresholds of dSmP ≤ 0.02 s and dSmP ≤ 0.01 L L\ns are used to ensure the 50% source overlap, respectively. We have added more details in Methods – Repeating Earthquake Identification (see lines 487-500).\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nFigure R6. Normalized distance (i.e., inter-event distance between repeater-candidate centroids divided by the sum of the rupture radii of the two events) as a function of dSmP over range of magnitudes. The stress drop is assumed to be 3.0 MPa in (a) and 19.1 MPa in (b), with Vp = 4 km/s.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nResponses to the comments by Reviewer 2\nQ 2-1 The manuscript written by Peng, Chen, and Burgmann describes their observations of the interaction of repeating earthquake-derived aseismic slip and seismic swarms leading to a large M 7.3 earthquake in eastern Taiwan. The overall presentation is clear and easy to understand. The data quality and presentations are highly reliable, providing a robust observation for the precursory slow slip before a large earthquake. This manuscript is significant for studying precursory slip leading to a large earthquake.\nBased on the observations, this manuscript discussed the possible controlling factors for earthquake swarms' diffusivity and fluid's origin. However, I do have one suggestion for the manuscript. Considering the 24-year catalog of repeating earthquakes, swarms, and regular earthquakes, at least five episodes of moderate to large earthquakes occurred in eastern Taiwan. With this dataset, a few parameters can be derived relating to their locations and timing: the aseismic slip at depth calculated from the repeating earthquake sequences, earthquake swarm diffusivity, and the occurrences of moderate to large earthquakes. The authors can derive an empirical relationship among these parameters and see how they relate. This relationship or formula may explain how the precursory slip leads to a large earthquake.\nResponse 2-1\nThe reviewer is interested in knowing if there exists a systematic spatiotemporal correlation between RES accelerations, swarms, and large earthquakes that can be expressed as an empirical formula. As shown in Figure R7a (also the current Figure S13), the separation distance between the swarms and the aseismic slip acceleration (defined by the peak in slip rate history in Figure 4) ranges from 10 to 50 km, corresponding to separation times of < 0.4 year; A spatial closeness of RES and swarms does not always produce a shortened interaction time, suggesting the driving mechanisms of RES and swarms are not necessarily the same.\nFigures R7b-c, on the other hand, demonstrate the correlation between M≥6 events and the inferred swarm/RES slip rate. We found that most of the earthquake swarms\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\noccurred at least 39 km from the M≥6 events. The swarm sequence with the shortest distance from an M≥6 event (27 km) is also characterized by the longest swarm duration and a short separation time, leading to an apparent positive correlation between the separation distance (Y) and time (X) (Figure R7b) as (1.425±1.486)X + (36.07±7.42).. Compared with swarm sequences, the acceleration of aseismic slip rate from RES reveals a much closer spatiotemporal association with M≥6 events, as revealed by their relatively short separation distances (5- 40 km) and times (< 0.35 yr) in Figure R7c. However, no clear scaling in the correlation is observed for RES’s spatiotemporal closeness from the seven M≥6 events. These results indicate that while swarms and M ≥ 6 events exhibit possible triggering relationships, RES accelerations appear more loosely coupled with subsequent large earthquakes.\nHowever, further statistical examination is needed to strengthen these correlation and interpretations. We have added this description to the Discussion, please see lines 198-210 for details.\nFigure R7 (new Figure S13). Spatiotemporal correlation between (a) swarms and aseismic slip rate inferred from RES, color coded by the swarm duration. The number in the circle represents the ID of the swarm sequence. (b) M≥6 events and swarms, color coded by the swarm duration. The number in the circle represents the ID of the swarm sequence. (c) M≥6 events and aseismic slip rate inferred from RES, color coded by the RES slip rate. The number in the circle represent the peak slip rate, while\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nthe number in the square represents the ID of M6 events. Note the horizontal axis in (a-c) indicates the absolute value of the time difference.\nResponses to the comments by Reviewer 3\nQ 3-1 This paper attempts to gain insight into aseismic processes related to the occurrence of Magnitude 6-7 class earthquakes by extracting repeating earthquakes and earthquake swarm activities from long-term seismic observation data. The authors claim to 'observe an acceleration of aseismic slip rate that occurred approximately three years preceding the April 3, 2024 Mw 7.3 Hualien earthquake'. While the proposed correlation between major earthquakes and aseismic slip is potentially very interesting, the supporting evidence presented in the current form falls short of being conclusive.\nThe main text is written based on the assumption that both repeating earthquake sequences (RES) and earthquake swarms are driven by aseismic slip (or slow slip), manifesting as seismic activity as RES and earthquake swarms. However, in general, the relationship between earthquake swarms and aseismic slip is not obvious and should be supported by additional observational data, such as geodetic measurements. Such data is not presented in the paper. Moreover, given the mention of the involvement of fluids in earthquake swarms, it is conceivable that fluid migration alone could explain earthquake swarm activity. The necessity of invoking aseismic slip in addition to fluid migration has not been adequately justified.\nResponse 3-1\nIn the revised manuscript, we have now focused on the period of accelerated aseismic slip from RESs and swarm activity (2019 to 2024) to explore both seismic and geodetic observations during the years preceding the 2024 Hualien earthquake. By analyzing Global Navigation Satellite System (GNSS) data from nearby continuously operating stations, we found that crustal deformation of 10-25 mm was observed during the swarm episodes in 2019 and 2021, suggesting the occurrence of aseismic slip during the swarm sequences, as shown above in Figures R4-5 (new Figure S14 and Figure 5). We have added a new section “4. Geodetic observation during the swarm sequences” in the revised manuscript (lines 259-333). Based on the clear migration\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\npatterns revealed by the swarm sequences, a combination of two driving mechanisms (propagation of fluids and slow slip) is proposed for the swarm activity in study area.\nQ 3-2 The discussion of hypocenter migration during swarm episodes is dependent on the definition and extraction methods of earthquake swarms. While Figure 4c plots only the extracted swarm activities, it is difficult for readers to judge the validity of the extraction within the overall seismic activity. For example, it would be necessary to show a magnitude vs. time plot similar to Fig. S2b within the rectangular area outlined in gray in Fig. 4a.\nResponse 3-2\nThe reviewer requested the visual inspection of background seismicity vs. swarm sequences during the 24-year study period, in order to judge the validity of seismicity declustering. Following this request, below we plot the magnitude and occurrence rate of the background seismicity against time (Figure R8). Different from the mainshock-aftershock sequence characterized by M6+ mainshock with accelerated occurrence rate, swarm sequences appear to coincide with the high occurrence rate with maximum magnitude in a range of M5.1~5.9 (except for the 2024 swarm on April 22). However, such visual inspection does not provide enough detail of the swarm extraction scheme used in this study. We think the reviewer’s concern can be better answered by showing the density rate (μ0) before and after applying individual declustering algorithms instead of plotting magnitude vs. time plot. We have added a new figure (Figure R9, same with Figure S25) and the corresponding description in Method “Earthquake swarm identification”. Please see lines 524-556 for more details.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nFigure R8. Magnitude and daily occurrence rate of seismicity as a function of time. Swarm events are colored with the same scheme in Figure 4.\nFigure R9 (new Figure S25). (a) Time series of cumulative background seismicity (ω ) 0\nfrom alternative declustering algorithms and cumulative earthquake numbers of raw data. Dashed vertical lines donate the occurrence times of M≥6.0 events. Black line indicates the raw data. Blue, red, and yellow lines indicate the declustered background seismicity using ETAS, Re85, and NN13, respectively. (b) Number of identified clusters using different minimum μ .\n0\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\n(c) Number of identified clusters with maximum magnitude ≥ 6 corresponding to minimum μ .\nResponse 3-3\nWe agree with the reviewer that the evaluation of migration of the 2021 swarm sequence can be improved by considering different representations of spatiotemporal plots that consider migration directions. The revised depictions lead to the following improvements: (1) They better describe the bilateral migration in the 2021 swarm sequence and (2) they clarify whether projecting the events onto movement direction (along-dip) provides better understanding of swarm migration and lend support to diffusion model.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nNow we have revised the figure by providing an along-strike cross-section with a better color scheme, which reveals a clear bilateral migration. This is shown in Figure R10 (also new Fig. 3) below. With the first event denoted by an open circle, this figure illustrates the upward and downward propagation directions along a similar dip-angle. Notably, 664 out of 966 swarm events (68.7 %) migrated upward, while the rest migrated downward.\nFigure R10. (a) Spatial distribution of the swarm sequence from April 18 to August 18, 2021. The circles indicate the swarm events, color-coded by the time elapsed since the beginning of the sequence. (b–c) Along-strike and across-strike distribution of swarm and RES events. (d) Magnitude-time plot of the swarm sequence. (e) Time-space evolution of the swarm sequence (see Supplemental Movie S2 for animated view of event sequence).The black solid line indicates the best-fit synthetic diffusion curve (D = 5.4 ± 1.1) determined by the majority of the 90th percentile distance points (pink squares) obtained each day, following 𝑟 =√4𝜋𝐷𝑡, where r is the distance (km) from the first event, t is the elapsed time (days), and D is the diffusivity (m2/s) with uncertainty denoted by one standard deviation (dashed lines).\nWhen we projected the hypocenters onto an inferred fault plane in a cross-section (Figure R11a), to calculate the along-dip distance by assuming propagation is primarily in dip direction, we found that the diffusivity (D) slightly changes with fault dip, in a range of 1.9 to 2.2 m2/s (Figure R11b). Considering a long-strike and along-depth distances reveal a similar range of D estimates from 1.2 to 2.7 m2/s. Compared to the\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nD estimated using just the hypocenter distance (~ 5 m2/s), the D measured using along-dip or along-strike distances are lower. This is mainly due to that comparing with the three-dimensional movement, the measure of two-dimensional expansion is likely underestimated. In addition, not every swarm sequence is characterized by clear migration patterns and thus it is difficult to find an optimal frame for projection axes, we prefer the standard approach to determine D using hypocenter separation. We have modified the original Figure 5 (now Figure 3) and provide more precise description of the migration pattern of the 2021 sequence. Please see lines 143-167.\nFigure R11. (a) Cross-section across the fault with varying dips of the fault plane. Orange, red, and blue lines indicate the assumed fault plane with dip-angle of 50°, 40°, and 35°. (b) Time-space plot illustrating all swarm events evolution using three fault-dip assumptions.\nQ 3-4 Furthermore, even if the diffusion model can be used for explaining the migration of the hypocenters in the extracted earthquake swarm activities, the diffusion model fitting shown in Figures 5c and S5-S13 is performed with respect to the distance from the first EQ of an extracted earthquake swarm cluster. Therefore, if a different earthquake within the same earthquake swarm activity had been chosen as the initial EQ, the spatial location of the reference for \"distance\" could change significantly, potentially leading to a large difference in the estimated value of D. From this perspective, as mentioned above, it is necessary to consider the examination along a distance axis in one direction within the 3D space.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nResponse 3-4\nWe fully agree that estimates of hydraulic diffusivity from distance-time plots of seismicity can be biased by choices of the origin and direction of propagation, and this is an issue with many published studies. To make our analysis more rigorous, a sensitivity test of D is now conducted by changing the choice of the initial events to provide a mean value of D with standard deviation. Here, instead of removing a single event multiple times, we removed the 90th percentile distance point (pink squares in Figure R12 below) one by one at the early stage of a sequence until the beginning 20 points are all removed. In other words, the first time we removed the first point (90th percentile obtained using the first 20 distance points) to measure D , the second time 1\nwe removed the initial two points for D , and the third time D for removing three points, 2 3\nand finally, D for removing 20 points at the first stage of the sequence. Finally, we are 20\nable to obtain the mean D and standard deviation using the 20 measurements of D (D 1\nto D ). The results are summarized in Table R1 (also added as the new Table S2). We 20\nfound that the D estimates increase with standard deviation, indicating a lower reliability for large values of D. The more focused swarm sequences (e.g., Seqs. 8 and 9) are characterized by the largest number of events in a sequence and the smallest standard deviation among the analyzed sequences.\nNote that in the previous version, the best fit curve was selected to span the majority of the 90th percentile distance points obtained each day. To obtain a smoother spatiotemporal tend, we now determine the 90th percentile every 20 data points. We followed the procedure of Amezawa et al. (2021), to find the minimum RMS misfit between the observed and predicted triggering front using moving time bins containing 20 earthquakes, with 10 events overlapping between consecutive bins. We have modified the text in lines 152-167 to describe the revised method and standard deviation of D using sensitivity tests.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nFigure R12. The hypocentral distance (in all direction) as a function of time from the first event in the 2021 sequence.\nTable R1. Detailed information of swarm sequences in this study\nQ 3-5 Regarding the case of Figure 5c, it could be argued that the sequence of seismic activity can be divided into three approximately distinct phases: an initial period of approximately 20 days, followed by a period of about 10 days starting around the 50th day, and a final period of about 40 days commencing around the 70th day. Even within each of these phases, hypocenter migration can be\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nobserved. It is therefore worth considering whether the possibility of such migration within these sub-clusters should be further investigated.\nResponse 3-5\nWe agree with the reviewer that there appear to be several phases of migration behavior in the 2021 swarm sequence. As the scope of this study is to characterize the spatiotemporal correlation between different types of aseismic slip indicators; that is, repeaters, earthquake swarms, and surface deformation (we have added GNSS analysis in the revised manuscript), we tend to not focus on the detailed patterns of migration of the swarms. Note that the uncertainty of earthquake location in Taiwan is approximately 3.1 ± 2.7, 1.3 ± 1.6, and 4.6 ± 3.9 km (Wu et al., 2013) in longitude, latitude, and depth respectively, which challenges further exploring the fine structure in a swarm sequence. A careful earthquake relocation should be conducted for nation-wide earthquake swarms in the near future.\nQ 3-6\nminor points:The authors set a minimum magnitude threshold of 2.0 for the earthquake catalog. A more in-depth discussion is warranted regarding the implications of this setting (Mc) on the analysis of earthquake occurrence rates. It is recommended that this topic be explicitly addressed in the Methods section.\nResponse 3-6\nAs shown by Figure R13 below, the Mc from 2000 to 2024 in the study area is estimated to be 2.0. We have added the description of earthquake catalog in Method (please see lines 477-485).\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nFigure R13. Calculation of the completeness magnitude (Mc) and b-value of the Gutenberg-Richter power law for the earthquake catalog in the study area from 2000 to March 2024, using the maximum likelihood estimation method. Squares and circles indicate the incremental and cumulative frequency-magnitude distributions, respectively.\nQ 3-7 Figure 1:\nThe focal mechanisms (beach balls), along with the text of earthquake year and Mw, are overlapping with symbols such as open circles and brown diamonds, making these symbols very difficult to see. The focal mechanisms and related text should be moved to a position further out, where they do not overlap with other symbols\nResponse 3-7\nWe have modified the figure accordingly, please see the modified Figure1.\nQ 3-8 In Figs. 3 and 4, the \"along strike\" distance is shown, but it is unclear where this corresponds to within Figure 1a. A distance axis should be added to Figure 1a to clarify this.\nResponse 3-8\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nWe have modified the label of Figure 3b from along strike distance to along strike (LVF) distance. Please see modified version.\nQ 3-9 Since there are several different symbols used in the figure, a legend should be included to explain their meanings. The black rectangular regions and the square bracket-like shape to the right of the letters A' and B' are not explained in the figure caption.\nSimilarly, in Figure 1b and 1c, the symbols are difficult to see due to overlapping text. There is also overlapping text related to different earthquakes, which is meaningless.\nResponse 3-9\nWe apologize for the missing label and overlapping text in the figure. We have added a legend to Figure 1 and replace the text with the event ID (corresponding to Table S1). Please see the modified Figure 1.\nQ 3-10 Given that Table S1 lists major earthquakes, it would be helpful to add numbers to these events and only show these numbers in the figure, or to use some other method to distinguish between different earthquakes.\nResponse 3-10\nWe have added event IDs of the M6+ events in Table S1 and used the ID numbers in Figure 1.\nQ 3-11 L82-89, Figure 2: It is unclear which range of earthquakes within the CRF and LVF directions (generally north-south) are sampled for the plot. Please indicate in Figure 1 the specific area for which the earthquake density is calculated.\nResponse 3-11\nWe appreciate reviewer’s comment on clarification of Figure 2. This figure is removed from the manuscript because the previous section of “seismic episodes showing LVF-CRF interactions” is considered not necessary (see Q1-5).\nQ 3-12 The density is described as \"each 1 square kilometer area and within 20-day time,\" but the color scale only uses \"low\" and \"high\" without specific numerical values. Please provide specific numerical values for the color scale.\nResponse 3-12\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nWe have modified the color bar by adding “event counts” corresponding to the density. Please see the modified version of this figure, which is now Figure S1.\nQ 3-13 The labels for major earthquakes (year, Mw) overlap with the density colors, especially for earthquakes after 2022, making it difficult to see the density. Please improve the visibility by repositioning the labels.\nIn L89, \"vertical red bar\" is mentioned, but it is unclear which part is being referred to. Additionally, there are several vertical blue dotted lines in the figure, but there is no explanation for these in the caption.\nResponse 3-13\nWe have modified the figure and the caption accordingly.\nQ 3-14 To demonstrate the simultaneous activation of seismic activity on both sides of the fault, I would suggest that the authors add a plot that shows the time-dependent cumulative number of earthquakes within defined areas on either side of the fault in a map similar to Figure 1.\nResponse 3-14\nWe have added a plot of cumulative number as a function of time plot, please see the revised Figure S1.\nQ 3-15 In addition to the overall 25-year period, zoomed-in plots for each of the five specific periods of interest should be included in the supplementary information for both Figure 2a (density plot) and the suggested cumulative number of earthquakes plot.\nResponse 3-15\nThe original purpose of Figure 2 was to focus on the interaction of seismic events. However, as suggested by the first reviewer (Q1-5), this aspect is less relevant and thus removed from the manuscript. In the revised manuscript, zoomed-in plots were designed to focus on post-2019 period instead (new Figure 6).\nQ 3-16:L85:The abbreviation \"LV\" is used before its definition is given in L115.\nResponse 3-16:\nWe added the definition of the LV accordingly (line 49).\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nQ 3-17 L101-104, Figure S1: Although the authors discuss aftershock activity in the LVF and CRF regions, while the cumulative seismic moment they present is interesting, it would be more common and informative to also show the time series of the cumulative number of earthquakes and magnitude-time plots. Additionally, in Figure S1, only the long-term changes over 25 years from 2000 are plotted, but for a more detailed evaluation of aftershock activity, it would be more appropriate to show zoomed-in plots for shorter time periods.\nAs pointed out for other figures, the presence of symbols like A-F and labels indicating the occurrence year of major earthquakes can obscure the plotted quantity (cumulative seismic moment). In particular, the \"D\" symbol obscures a line, which is unacceptable.\nResponse 3-17\nThe previous Figure S1 is no longer needed in the revised manuscript, as we have decided to remove the interactive seismic episodes from the main text, as suggested by reviewer 1 (Q1-5).\nQ 3-18 Movie 1:\nI initially assumed that the color of the symbols represented the occurrence time of earthquakes. However, the presence of additional colors like gray and black makes it difficult to interpret the data. Furthermore, the limited colors for the circle symbols, primarily pink and blue, making it difficult to distinguish the differences in occurrence times. The overall purpose of this movie is unclear.\nMoreover, despite being a movie format, the lack of time or date labels for each frame makes it very difficult to understand what the author wants to show. I think it would be more clear if the time scale in the lower left corner had an arrow indicating the corresponding time, and the arrow moved to the right as time progresses.\nResponse 3-18\nWe agree that the current movie can be improved. Now we have modified the movie to illustrated the sequential ruptures from 2021 to May 2024, including the main active zone of the earthquake swarm (west-dipping CRF), main seismogenic zone on the NLV (2022 Mw 6.1 and its associated aftershocks, west-dipping CRF), and early aftershocks of the 2024 Mw7.3 earthquake (east-dipping LVF).\nQ 3-19 Figure 3a:\nAs I comment on Figure 1, the text labels indicating year and Mw of\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nearthquakes overlaps with other symbols, making them very difficult to see.\nThe figure caption mentions \"shaded and open polygons,\" however, it is not clear which specific features of the figure these terms refer to.\nResponse 3-19\nWe have modified the figure by replacing the time of events by their event ID and removed the original Figure 3a.\nQ 3-20 Figure 3b:I believe that overlaying a color map of earthquake density with the RES slip rate would compromise the clarity of the figure. I recommend presenting these two datasets separately.\nFurthermore, the vertical axis of the RES slip rate plot should be clearly labeled, particularly for the creeping CRF slip rate plotted at the top. The current plot makes it difficult to identify the zero level (0 cm/day), and the data appears to extend into negative values. Additionally, the labels indicating major earthquakes overlap with the data lines, especially after 2022, making it difficult to see. These issues should be addressed before finalizing the figure.\nResponse 3-20\nWe have modified the figure accordingly, please see the revised version of Figure S1.\nQ 3-21 The color scale in this figure is described as representing \"event density\" in the figure caption. However, given the presence of the RES slip rate within the same figure, this could be interpreted as referring to the event density of repeating earthquakes. To avoid any ambiguity, please provide a more specific description that clearly distinguishes between these two quantities.\nResponse 3-21\nWe apologize for the confusing label. We have modified the figure by separating the RES slip rate from the density plot.\nQ 3-22 Figure S2b: The cumulative number of EQs and the magnitude-time plot overlap, and the labels indicating the years and Mw of major earthquakes, as well as the colored sections for different time periods, further obscure the background magnitude-time plot. As mentioned in other figures, please improve visibility by plotting these elements separately.\nResponse 3-22\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nWe have modified the plot accordingly. Please see the revised Figure S2.\nQ 3-23 L156-159:\nWhile the author states a correlation between the aseismic slip rate and RES occurrence rate around 2003-2004, a direct comparison of Figure 3b and Figure S3b is required to verify this claim. However, both figures are plotted over very long periods, making it difficult to compare the changes in the 2-3 year period mentioned here. I would recommend creating a new figure with Figure 3b and Figure S3b plotted vertically side-by-side, sharing a common x-axis (time axis) and zooming in on the 2-3 year period of interest.\nResponse 3-23\nThank you for your comments. The comparison between the aseismic slip rate and RES occurrence rate is shown by the figure below (new Figure 2). In the southern LV area, an increasing number of repeating events (vertical bars in the upper panel) occurred immediately after the 6.4 mainshock (15 repeating events in one week), resulting in a sudden change of RES occurrence rate (black line in the bottom panel). The statement of “The highest aseismic slip rate peak occurred several months after the 2003 Mw6.8 Chihshang event, coinciding with an acceleration in the RES occurrence rate” is dependent on the averaging window adopted for slip rate estimates. In order to show long periods characteristics, the averaging window is selected as 0.5 year in Figure R14 below (new Figure 2), where the RES slip rate was calculated as the average value within a 3-month window before and after each plotted point. A shorter moving window provides higher temporal resolution but may be less meaningful if it is shorter than the interevent time in RES catalog. A longer moving window is more statistically meaningful but may produce time shift of an actual signal. In this study, the mean and maximum interevent times are 15 and 127 days, therefore the calculation of slip rate is less reliable using time window shorter than 4 months. Particularly, the near-zero slip rate could be a result of no data point. We have discussed the contribution of time window selection in the temporal characteristics of RES-derived slip rates in lines 103-109.\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nFigure R14 (new Figure 2). Magnitude distribution (top), occurrence rate and slip rate (bottom)of RESs since 2000 for (a) northern LV on the CRF and (b) southern LV on the LVF. Cumulative RESs events and the derived slip rate are denoted by black and dark red lines, respectively. Vertical red dashed line indicates the occurrence time of Mw ≥ 6 events with the ID number corresponding to Table S1. Black arrow and black dashed line mark the significant change of the RESs occurrence rate. The occurrence rate averaged over each sub-period is indicated by the number.\nQ 3-24 Figure S3:\nThe M-T plot and the cumulative RE number line overlap, making it very difficult to see. To improve clarity, I suggest dividing the figure into two panels, sharing a common x-axis.\nWhile the figure caption mentions \"The green, red, blue sections,\" the use of color alone to represent different time periods may be difficult to distinguish. I think it would be better to use a shape like a square bracket (]) at the top of the top panel to indicate the time periods.\nAdditionally, please provide a detailed explanation in the figure caption for all elements, including arrows, red dashed lines, and numerical values.\nResponse 3-24\nThe original Figure S3 illustrated the temporal distribution of the occurrence rate and magnitude of RESs. We have modified the figure accordingly based on the reviewer’s suggestion. Now the revised figure is shown in the new Figure 2.\nQ 3-25 L199-200:The authors excluded swarm sequence #10 from their discussion\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nbecause they have not updated the RES catalog.\nHowever, I think it would be better to update the catalog and include the swarm sequence #10 in the analysis. The current approach seems rather incomplete.\nResponse 3-25\nWe do not yet have the waveform data recorded at the same stations used in the original RES identification. Thus, within the limited time frame for the revision, we are not able to provide this update. We suggest that this does not really affect the strength of our results and conclusions.\nQ 3-26 L200-201:\nThe authors states that \"6 out of 9 swarm sequences were preceded by the increasing aseismic slip rate,\" but Figure 4b shows that only #1, #2, #3, and #8 exhibit this pattern. For sequences #7 and #9, the aseismic slip rate and the swarm onset appear to occur simultaneously.\nAdditionally, the confidence interval of the RES slip rate is also unclear.\nFurthermore, the selection of swarms is dependent on the definition used in this study. The robustness of such arguments should be examined if the definition is changed.\nResponse 3-26\nWe appreciate this suggestion. To address the question of 'how aseismic slip is associated with the swarm activity', we have added a new section “Spatiotemporal association between swarms and repeating earthquakes” where we performed a sensitivity test (Figure S12) by computing the averaged RES slip rate before (d ) and pre\nafter (d ) each swarm onset within varying time windows (1-12 months). When the post\nslip rate is averaged over the time interval of 3 to 6 months, we found a consistent pattern with 6 out of 9 swarm sequences preceded by an increasing aseismic slip rate. Seqs. 2 and 9 exhibit the highest spatial overlap with RESs (light green squares in Figures S4 and S11), corresponding to the strongest temporal coincidence, as indicated by the largest RES slip rate change across the swarm sequence (Figure S12). Please see lines 186-197 for more details.\nQ 3-27 Figure 4:\n(a) Similar to the other figures, the text overlaps with the symbols.\n(c) The overlay of swarm locations (along strike distance) and the RES slip rate results in portions of the black line being obscured. To improve clarity, the plotted swarm locations should be adjusted vertically to avoid overlapping with the black line.\nResponse 3-27\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nWe have modified the plot accordingly. Please see the revised Figure 4c.\nQ 3-28 Figure 4: It should be clearly stated in the caption whether the plotted RES slip rate (black line) is identical to one of the two RES slip rates plotted in Figure 3b, or if it represents a different value. The current state of the figure is likely to confuse readers.\nResponse 3-28\nThe RES data used for slip rate estimate in Figure 4c is denoted by grey box in Figure 4a, lightly smaller than the one for old Figure 3a (now Figure S1). We have clarified the difference in the caption of Figure 4c.\nQ 3-29 The diameters of the open circles vary, but there is no explanation for what these differences represent (similarly in Figure 1).\nResponse 3-29\nWe have added the legend for the size of circles, please see the revised version of Figure 1, Figure 3 and Figure 4.\nQ 3-30 Movie 2:\nAs pointed out in Movie 1, it would be better to show the time of each frame in the movie.\nSimilarly, the color scheme for the symbols should be improved as I pointed out in Movie 1. The sudden appearance of numerous points at the beginning of the movie creates a strange impression as it is unclear how simultaneous they are because the time interval between each frame is unclear.\nResponse 3-30\nWe have modified the movie accordingly.\nQ 3-31 L219:\nThe text mentions \"low Vp and Vs,\" but this data is not presented in any of the figures. At least a citation should be provided.\nResponse 3-31\nWe added the citation, accordingly.\nQ 3-32 L277: \"located to the south of the earthquake swarms, opposite to their propagation direction\".\nPrior to this sentence, the text has described the vertical migration of hypocenters from depth to a shallower part. The abrupt introduction of the term \"south\" seems to imply a migration in north-south direction. This ambiguity makes\nWei Peng1, Kate Huihsuan Chen1, Roland Bürgmann2, Ya-Ju Hsu3, Yan-Hong Chen3\nit difficult to determine the exact direction of migration.\nResponse 3-32 By adding the along-strike cross-section in old Figure 5, the northward propagation is now better illustrated. Please see the modified text associated with new Figure 3.\nQ 3-32 Figure 6a:\nAlthough four quantities are plotted in this figure, it is very difficult to understand due to the unclear vertical axes. To improve clarity, I suggest dividing the data into four separate subplots, each with a clear vertical axis and a shared horizontal time axis.\nResponse 3-32 We have modified the figure accordingly. Please see the revised Figure 6a to 6e.\nQ 3-32 Figure S4:\nSame comments as for Figure 1.\nResponse 3-32\nWe have modified the figure accordingly.\nQ 3-32 Figures S5-S13: Same comments as for Figure 5.\nResponse 3-32\nWe have modified the figures accordingly.\n32\nWe appreciate the reviewers’ construc2ve sugges2ons to improve this manuscript.\nIn this revised version, we have strengthened the presenta2on of the accumula2ng\naseismic and seismic phenomena since 2023 by incorpora2ng detailed analyses of\nrepea2ng earthquake moment rates, a- and b-values, and declustered seismicity\nrates (modified Figure 6). Addi2onally, we have added a new Table S5 and example\nin Figure S26 to describe the correc2on of aMerslip following individual M6+ events.\nThe point-by-point responses to each comment are shown below.\nReviewer #4 (Remarks to the Author):\nQ 4-1 This paper investigates the process leading to the Mw7.3 earthquake in the double-\nvergence suture zone in eastern Taiwan based on repeating earthquake sequences,\nearthquake swarms, crustal deformation, and static stress change analyses. The\nresults suggest that aseismic slip events and transient fluid flow potentially promoted the\n2024 mainshock. The revised manuscript has responded sincerely to the reviewer’s\ncomments, and the reliability of the data analysis appears to have been significantly\nimproved. The contribution of aseismic slip and fluid flow to each earthquake swarm and\nrepeating earthquake activity is well described in the main text. The finding that stress\nincreases due to these phenomena may have promoted the occurrence of large earthquakes\nis of great significance. However, this manuscript does not provide clear evidence for the\nacceleration of seismic and aseismic slip rates about one year prior to the 2024 M7.3\nearthquake (Stage 3: Lines 306-315). Shallow earthquake swarms do not necessarily\nappear to occur at the fault plane of the CRF. It is not clear from the illustrated repeating\nsequences and GNSS data whether the CRF was creeping at a relatively fast rate. I\nrecommend that the process for this period be described more clearly.\nAns:\nWe thank the reviewer for pointing out the need to better clarify the evidence for acceleration\nin seismic and aseismic slip rates approximately one year prior to the 2024 Mw7.3 Hualien\nmainshock (Stage 3). To address this, we now present three lines of evidence that together\nsupport a phase of evolving deformation on the Central Range Fault (CRF):\n(1) Increase of RES activity\nAs shown in Figure R1, we highlight the observed increase in the seismic moment rate of\nrepeating earthquake sequences (RESs) across 2023, which rose from approximately 1.3–1.4\n× 10¹⁵ to 2.2 × 10¹⁵ N-m/day, while the number of RES events remained relatively stable (~45\nevents/year). The observed increase in the seismic moment of RES one year prior to the 2024\nMw7.3 coincides with a minor but consistent increase in RES-derived slip rate in the long-term\nmoving averages (black line in Figures 6b), although it may not stand out in short-period\nwindows (green line in Figures 6b). Such ongoing minor acceleration of aseismic slip may not\nbe fully captured by GNSS data (Figure 6e) due to the depth of the RESs (~15–25 km) and\nsparse geodetic coverage in the mountainous hanging wall.\nThese observations support a scenario where aseismic slip and associated stress loading of\nasperities intensified across 2023, leading to a rise in RESs seismic moment without a\nproportional increase in their number. This progressive fault weakening, potentially aided by\ndeep fluid migration and undetected deep slow slip, may have played a crucial role in\npreparing the CRF for dynamic rupture, ultimately contributing to the nucleation of the April\n3, 2024 Mw7.3 earthquake.\nFigure R1 (same with the new Figure 6b-d). The temporal distribution of RES’s magnitude\n(vertical grey bar in (a)), slip rate (black line in (a)), cumulative Mo (black line in (b)), and\ncumulative number (black line in (c)).\n(2) Gradual increase in background seismicity\nTo examine the evolution of background (declustered) seismicity, we analyzed both shallow\n(≤15 km) and deep (>15 km) earthquakes (Figure R2). We found that shallow background\nseismicity began a clear upward trend approximately one year before the mainshock, while\ndeeper seismicity remained relatively stable. This may suggest that a progressive stress\naccumulation or fault weakening were concentrating at the shallow depths.\nFigure R2 (same with Figure 6i). Declustered seismicity at shallower (dark blue) and deeper\n(light blue) depths.\n(3) Increase in a- and b-values\nWe also computed the time-varying a- and b-values using the original background seismicity.\nAs shown in Figure R3, both a- and b-values began increasing in the shallow seismicity\nfollowing the cascade of M6+ events since 2023. The rising a-value reflects an overall increase\nin seismic activity, while the increasing b-value suggests enhanced participation of smaller\nevents. This is consistent with distributed micro-failure and increasing heterogeneity or fluid\neffects. In contrast, deep seismicity exhibited no notable changes. Such depth-dependency is\nconsistent with the declustered seismicity in Figure R2.\nTo reconcile the apparent increase in deep aseismic slip and RES moment rate with the\nshallow increase in declustered seismicity rate and in a- and b-values, we propose a depth-\ndependent fault evolution model. At greater depths (>15 km), the CRF is undergoing\nprogressive aseismic slip, potentially accompanied by upward fluid migration, may contribute\nto a weakening of the upper fault zone. The increase in shallow seismicity, coupled with higher\na- and b-values, suggests a brittle response from a more heterogeneous and potentially fluid-\ninfluenced fault environment. These observations suggest that deep slow slip and shallow\nbrittle failure operated as coupled processes along the CRF, representing an intermediate\nstage in the preparatory phase leading to large earthquake nucleation. We have modified\nFigure 6 and changed the texts for Stage 3 accordingly, please see lines 311-348 for details.\nLastly, to address whether the swarm sequences occurred on the fault plane of the CRF, we\nrefer to Figure S2 for clarificafon. The projecfon of all idenffied swarm events shows that the\nmajority of events align with and overlap the west-dipping structure, consistent with the\ngeometry of the CRF.\nOther minor corrections in the manuscript are listed below.\nQ 4-2 Figure 3: The figure is duplicated. Please confirm and delete the unnecessary figure.\nAns: We apologize for the mistake; we removed the unnecessary figure accordingly.\nQ 4-3 Figure S1 caption: (b) … for the southern LV area -> for the northern LV area\nAns: We have modified the capfon accordingly.\nQ 4-4 Figure S15: The arrows for predicted displacement and observed displacement are\ndifficult to distinguish. Please change to a more visible color.\nAns: We have modified the color accordingly, please see the new Figure S15 (Figure R3 below).\nFigure R3 (same wit Figure S15). The opfmal forward models during the (a) 2019 and (b) 2021\nswarm sequences. Orange lines indicate the top of the fault. Blue vectors indicate the imposed\nslip. Red and black vectors indicate the predicted and observed displacement, respecfvely.\nQ 4-5 Figure S22: The units in the figure should be changed to SI units (kPa instead of bar).\nAns: We have modified the unit to kPa accordingly. Please see the modified Figure S22.\nQ 4-6 Several errors were found in the occurrence year and magnitude of the target\nearthquakes and in the figures cited, which should be corrected.\nAns: We appreciate the reviewer’s careful reading and aoenfon to detail. In response, we have\nthoroughly reviewed the enfre manuscript, figures, and supplementary materials to idenffy\nand correct all typos and errors, including inaccurate earthquake years, magnitudes, and figure\ncitafons.\nQ 4-7 Line 61: Movie S1 -> Movie S2\nAns: We have modified the typo accordingly.\nQ 4-8 Line 66: Figure 1a and 1b -> Figure 1b and 1d\nAns: We have modified the typos accordingly.\nQ 4-9 Line 140: Figure 1a, Figure S2a-c -> Figure 1b-c, Figure S2\nAns: We have modified the errors accordingly.\nQ 4-10 Line 168: table S3 -> Table S2\nAns: We have modified the typo accordingly.\nQ 4-11 Line 227: M6.6 -> M6.1\nAns: We have modified the typo accordingly.\nQ 4-12 Line 242: Movie S2 -> Movie S1\nAns: We have modified the typo accordingly.\nQ 4-13 Lines 286, Figure 6 legends and caption: “< 15km, ≥ 15km” or “≤ 15 km, > 15km”\nAns: We have modified the typos accordingly.\nQ 4-14 Line 324: Figure 6b -> Figure 6f\nAns: We have modified the typo accordingly.\nQ 4-15 Line 346: April 22, 2023 -> March 22, 2022\nAns: We have modified the typo accordingly.\nQ 4-16 Line 360: Figure 22 -> Figure S22\nAns: We have modified the typo accordingly.\nQ 4-17 Figure 5: (a) 2021 Mw5.8 -> 2021 Mw5.9\nAns: We have modified the typo accordingly.\nQ 4-18 Line 871 (Figure 5 caption): Figure 1a -> Figure 1b\nAns: We have modified the typo accordingly.\nQ 4-19 Line 891 (Figure 6 caption): gray box in Fig. 3a -> Fig.4a\nAns: We have modified the error accordingly.\nQ 4-20 Figure S18 caption: 2023 -> 2022\nAns: We have modified the typo accordingly.\nQ 4-21 Section 4 had different descriptions in the merged file and track change file. Please\nconfirm which is the final version.\nAns: We apologized for the mistake; the merged file was the final version. In the new\nsubmioed files we will avoid the same mistake.\nReviewer #5 (Remarks to the Author):\nQ 5-1 Are “Repeaters” and “RSEs” referring to the same events? If so, consider using a single\nterm consistently throughout the manuscript to avoid confusion.\nAns: We modified the wording from repeaters to RESs throughout the manuscript for\nconsistency, thank you for the suggesfon.\nQ 5-2 Figure 3. Please provide clearer details on the color-coding scheme used in Figure3. The\ncurrent figure does not clearly explain what each color represents.\nThe version of Figure 3 in the response letter appears different from the one in the manuscript.\nAdditionally, the merged manuscript includes two sets of (a), (b), (c) subfigures under Figure\n3. Please confirm and correct as necessary.\nAns: We have modified the figure accordingly, please see the new Figure 3.\nQ 5-3 Response 3-22: The authors mention modifying the magnitude-time (M-T) plot in Figure\nS2. However, the changes are not clear and do not seem to fully address the reviewer’s\noriginal concerns.\nAns: We apologized for the mistake in the rebuttal. The original Figure S2 (Figure R4 below)\nwas removed from the manuscript. It is no longer used in the current manuscript.\nFigure R4. (a) Distribufon of earthquake acfvity since 2020. Three groups of avershocks\nwithin one month of the M6+ mainshocks are indicated by yellow, orange, and purple circles\nfor the 2022 M6.7, 2022 M6.1, and 2022 M7 events, respecfvely. (b) Number of earthquakes\nin the selected area denoted by the dashed square in (a). The avershocks within one month\nof the three mainshocks are indicated by the same colors as in (a). The green line indicates\nthe cumulafve number of earthquakes. The dashed lines highlight the dramafc change in\nearthquake occurrence rate following the March 22 2022 M6.7 event. Note that the\nearthquakes plooed here are from the original catalog without declustering.\nQ 5-4 Response 3-23 and Figure 2: The black arrow between 2011 and 2012 in the bottom\npanel of subfigure (b) does not appear to coincide with a notable increase in RSE occurrence\nrate. What is your criterion for defining a “significant” change in RES occurrence rate? In the\nmanuscript or caption, can authors provide more details on the sub-periods used to calculate\nRES occurrence rate?\nAns: We appreciate the reviewer’s observation. The criteria for defining a “significant” change\nin RES occurrence rate were not clearly explained in the original manuscript or figure caption.\nTo clarify, the three phases shown in both panels (a) and (b) of Figure 2 are defined based on\na relative change of at least 20% in the RES occurrence rate compared to the previous time\nperiod. For example, in panel (a), the transition from 41.7 to 33.3 events/year represents a\n~20% decrease, and the subsequent increase to 45.8 events/year corresponds to a 37% rise.\nIn panel (b), the change is more pronounced: from 13.7 to 41.6 events/year (~203% increase),\nfollowed by a decrease to 30.5 events/year (~27% decrease). To clarify, we have added the\ncorresponding description in the caption of Figure 2.\nQ 5-5 Movie S1. While you indicate that circle colors represent days since April 18, 2021,\nconsider adding a color bar legend to clarify this temporal coding across frames. It is also\nunclear whether RES locations are marked in the animation. Please clarify how they are\nrepresented and revise the caption accordingly for better interpretability.\nAns: Thanks for the comment. The swarm events were color coded by the days aver the first\nevent (April 18, 2021), as plooed on the lev penal. In the revised version, we have enlarged\nthe legend on the lev penal. We have also added repeafng earthquakes in Movie S1, to show\ntheir spafotemporal associafon with swarm events.\nQ 5-6 Figure S4. The figure caption is confusing. “As for Figure S4” seems misplaced, should\nthis refer to Figure S3? Please review and ensure clarity across the following similar figures.\nAns: We apologize for the mistake, we have modified the capfon from “as for Figure S4” to\n“as for Figure S3”, accordingly. Please see the new version of Figures S4-S11.\nQ 5-7 Line 132: Please explain “10 km of each other”.\nAns: The “10 km” criterion is based on our previous study (Peng et al., 2021), where we\nsystemafcally examined the spafotemporal relafonship between M ≥ 3 swarm sequences\nand M ≥ 6 earthquakes in Taiwan. This threshold was adopted to assess whether short\ntemporal separafons between swarms and larger earthquakes are also associated with small\nspafal separafons. Our analysis revealed that a significantly higher proporfon of swarm\nevents tend to occur following M ≥ 6 earthquakes when the separafon distance is less than\n10 km. This finding suggests that swarm acfvity is more likely to be influenced by nearby large\nearthquakes when occurring within this proximity, as illustrated in Figure R5 below. We have\nadded descripfon to explain the “10 km” with citafon, please see lines 130-132.\n[REDACTED]\nFigure R5 (same with Figure 7a in Peng et al. (2021)). Relafve fming of swarm events in\nrelafon to M ≥ 6 earthquakes across varying spafal separafons. The distribufon of fme\ndifferences between 4,726 swarm events and 59 M ≥ 6 earthquakes. Each fme difference is\ncalculated as the occurrence fme of the swarm event minus that of the M ≥ 6 event. The cool-\ncolored lines represent expected distribufons from 50 synthefc M ≥ 6 catalogs simulafng\nrandom temporal associafons. This figure illustrates how the likelihood of swarm acfvity\nincreases aver large earthquakes, especially when occurring within shorter spafal distances.\nQ 5-8 Line 53: The phrase of “there systems” is unclear. Please clarify.\nAns: We have modified “there systems” to “these fault systems” for beoer refers to the \"two\nopposing faults parallel, head-to-head fault structures\" menfoned in the previous sentence.\nPlease see the line 53.\nQ 5-9 Line 105: The statement “We also test sensitivity to different time windows” is vague.\nPlease clarify what variable or outcome the sensitivity refers to.\nAns: Thanks for the suggesfon, we have changed the phases as “We systemafcally tested the\nsensifvity of the computed slip rate” please see the lines 103-104.\nQ 5-10 Lines 135-137: This sentence is unclear, especially the later part. Please rewrite this\nsentence.\nAns: We have modified the sentence to “Swarm sequences were idenffied from events with\nM ≥ 2.0, consistent with the esfmated magnitude of completeness (Mc) for this region“.\nPlease see the new Lines 134-135.\nQ 5-11 Line 140: “There northern swarms” should likely be “The northern swarms”.\nAns: We have fixed the typo accordingly.\nQ 5-12 Line 186: “To test the robustness of this relationship” – which specific relationship are\nyou referring to? Please be explicit.\nAns: We have modified the phrase by adding “between RESs accelerafon, swarms, and M≥6\nevents”.\nQ 5-13 Lines 216-217: Are postseismic deformation signals evident in the GNSS data? The\nMethods section describes modeling postseismic relaxation using an exponential function-\nwhy was this specific function chosen, and why is postseismic effect not discussed in the\nResults?\nAns: We thank the reviewer for highlighting the importance of clarifying the treatment of\npostseismic deformation in our analysis. Although we did not explicitly discuss postseismic\ndeformation in the Results (since our primary focus was on aseismic slip signals during swarm\nsequences), we have carefully accounted for postseismic effects during GNSS data processing,\nwhich is now described in the Methods section.\nWe applied a station-specific least-squares regression approach to model and correct for\npostseismic relaxation at sites showing clear coseismic offsets, including those affected by the\nexponential decay function of the form\n∆(𝑡) = 𝐴𝑒!\"/$ \n, where the parameters were calibrated through visual inspection to capture the observed\nrelaxation behavior. The example of fitting model versus observation is shown in Figure R6\nbelow. Please see the texts in lines 628-645 for details.\nFigure R6. Example of postseismic relaxafon in the north–south (N–S) component of GNSS\ndisplacemen. Red line represents the modeled trend considering long-term, seasonal trend,\nco-seismic offset and postseismic effect. Nearly verfcal red line (on the top of blue dashed\nline) indicates the fming of coseismic offsets corresponding to 2019 M6.1 event. The gradual\ndecay following the coseismic step reflects postseismic deformafon, which is modeled and\nremoved using an exponenfal funcfon during GNSS fme series processing.\nQ 5-14 Line 245: “die to ” -> “tie to”.\nAns: We have modified the world accordingly.\nQ 5-15 Line 230: Can authors provide further discussion on the potential slow slip or afterslip\nof the M5.9 event? How does this relate to subsequent large magnitude earthquakes?\nAns: We thank the reviewer for this constructive suggestion. As shown by Figure 5b-c,\ntransient deformation is evident during swarm activities. The 2019 swarm, which occurred\napproximately four months after the Mw 6.1 earthquake and ~30 km to its south, produced\n~10 mm of eastward and southward displacement at stations FLNM and DNFU. During the\n2021 swarm sequence particularly between Mw 5.9 and Mw 5.5 event, a gradual eastward\ndisplacements of ~ 10 mm were observed at stations YENL, TUNM, SOFN and up to ~25 mm\nat NDHU.\nWe interpret these displacements as primarily reflecting a combination of post-seismic and\naseismic deformation, rather than coseismic slip from the Mw 5.9 event at the depth of 15\nkm and its afterslip. This is supported by the concentration of seismicity between the two\nevents (M5.9 and M5.5) at depths around 10 km (Figure 3), which are unlikely to produce\nobservable surface displacements. Moreover, the GNSS displacement time series during this\ninterval can be well modeled by a post-seismic exponential decay function without requiring\na co-seismic step, as shown in Figure R7 and Figure S26. Please also see Method and the lines\n223-234 for details.\nRegarding the reviewer’s question about whether this postseismic deformation may relate to\nthe subsequent Mw 5.5 event, we acknowledge that addressing this would require a fully\nresolved fault model. However, as shown in Figure S15, the number and distribution of GNSS\nstations with detectable displacement are insufficient to robustly constrain such a model. We\nconsider this a promising direction for future work once denser geodetic coverage becomes\navailable.\nFigure R7 (same with Figure S26). Example of the GNSS fme series versus modeled\ndisplacement curves. (a) The east components of GNSS displacements at Three stafons (NDHU,\nTENL, SOFN), with black dots represenfng observed data and red lines denofng the fioed\ncurves, which include linear trends, seasonal signals, and postseismic exponenfal decay\nfuncfon aver M5.9 event. Verfcal lines mark the fming of the 2021 M5.9 and M5.5 swarm\nevents. (b) Mapview showing the epicenters and focal mechanisms of two targeted events\n(colored stars), GNSS stafon locafons (triangles), and earthquake swarm distribufon in 2021\n(black dots).\nQ 5-16 Line 263: Why authors show raw GNSS data instead of processed data?\nA: Sorry for misleading wording, we actually used the corrected GNSS data in old Figure 6e, as\ndescribed in the figure capfon.\nQ 5-17 In section 6. How do you define the significance of static Coulomb stress changes? Is\nany positive stress change considered significant? A stress change of 2-3 kPa is not large – how\ncan authors be confident that such a value reflects a considerable effect from aseismic slip?\nAns: At the end of Secfon 6, we menfoned that “the combined effect of both aseismic and\nseismic events likely plated a crucial role in the stress triggering of the 2024 M7.3 event. This\nis supported by the cumulafve stress increase of up to 30 kPa at the M7.3 rupture locafon,\ngenerated by aseismic episodes in 2021 along with three major ruptures in 2022”. We\napologize for not making it very clear. We have strengthened this paragraph by the following\ntexts below:\nInstead, the combined effect of both aseismic and seismic events likely played a crucial role in\nstress triggering of the 2024 Mw7.3 earthquake. This interpretafon is supported by the\ncumulafve Coulomb stress increase of up to 30 kPa at the eventual rupture locafon, resulfng\nfrom deep aseismic slip episodes in 2021 and three major M6+ ruptures in 2022 (Figure S22).\nNotably, a stress change on the order of 10 kPa is commonly considered sufficient to influence\nfault failure fming in tectonically loaded systems. Therefore, a cumulafve 30 kPa increase\nrepresents a substanfally elevated stress state, capable of significantly advancing the fming\nof a large earthquake. In addifon to fault weakening due to upward fluid migrafon, this level\nof stress accumulafon from both slow and fast slip along the CRF likely contributed to the\nenhanced seismicity and RES acfvity observed during Stage 3, ulfmately facilitafng\nnucleafon of the April 3, 2024 Mw7.3 mainshock. Please see lines 401-407.\nReviewer #4 (Remarks to the Author):\nI appreciate the author’s efforts to revise the manuscript. I think this manuscript will\nbe acceptable after some minor corrections have been made.\nLines 213-219: The same sentences are repeated. Please correct them.\nAns:\nWe have deleted the repeated sentences.\nLines 314-317: If the magnitude of a repeating earthquake in the same group was\nlarger than that of other earthquakes, it would be a significant change, but no such\nevidence is shown here. Although various factors should be considered, I think\n“significantly” is an overstatement in any case.\nAns:\nWe agree with the reviewer’s comment. We have replaced \"significantly\" by “notable”.\nLines 328-335: Please confirm whether the changes shown in this paragraph are\nresults obtained from earthquake information before the M7.3 mainshock. If the\ninformation includes data after the M7.3 mainshock, it should be described in\nStep 4. The 2024 swarm activity began on April 22, 19 days after the mainshock.\nTherefore, I am not sure it is appropriate to refer to Figure S17 in Line 329.\nAns:\nThis paragraph should have been removed because the anomalous a- and b-values\nobserved 20 days prior to the mainshock may not be statistically convincing. We had\ndecided to remove it from the text, but it was mistakenly retained in the revised\nversion (and Figure S17 was therefore not the figure used in the current\nsupplementary file). We sincerely apologize for this mistake and appreciate the\nreviewer’s attention to this issue.\nSeveral errors were found in the figures cited, which should be corrected.\nLine 109: Figure 1c -> Figure 1d\nAns:\nWe have modified from Figure 1c to Figure 1d, accordingly.\nFigure 3 caption: Please add an explanation for figure (d) and change “(c)” to “(e)”.\nAns:\nWe are not certain what the reviewer is referring to; however, after carefully reviewing\nthe caption, we did not identify any apparent issues. Nonetheless, we have added\nexplanations for the symbols to the caption for clarity.\nFigure 4 caption: Figure 1a -> Figure 1b\nAns:\nWe have modified the caption accordingly.\nFigure 6 caption: Fig. 3b -> Fig. 4a\nAns:\nWe have modified the caption accordingly.\nFigure 7 caption: Movie S2 -> Movie S1\nAns:\nWe have modified the caption accordingly.\nFigure S3 caption: 2020 -> 2000\nAns:\nWe have modified the caption accordingly."}, "author_response": ""}]

Decision letter