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 Mw7.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 Mw7.3 event.

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. In eastern Taiwan, RESs have been documented along the creeping segments of the LVF and CRF. We updated the previously established RESs catalog 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, 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, the slip for each repeating event can be obtained using the scaling relationship between the slip estimate (d) and seismic moment (Mo) as logd = - 1.21 + 0.11 logMo and logd = - 1.96 + 0.14 logMo 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.

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. 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 Mw6.8 Chihshang earthquake, coinciding with a sharp increase in RES activity. Following the Mw6.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. The most pronounced CRF aseismic slip acceleration followed the April 18, 2019 Mw6.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 a typical aftershock decay pattern. These swarms are generally associated with fluid intrusions and elevated pore pressure e.g., refs. , which can induce fault unclamping and promote slow slip. A previous study in Taiwan 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, 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 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 M2.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 Mw7.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. 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 , 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. The best-fit curve was selected to span the majority of 90 percentile distance points (pink squares in Fig. 3e) obtained every 20 earthquakes. 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 (90 percentile obtained using the first 20 events) to measure D, the second time we removed the initial two points for D, and the third time D for removing three points, and finally, D for removing the first 20 points at the earlier stage of the sequence. This yielded a distribution of diffusivity estimates (D to D), 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 Mw7.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 (d) and after (d) each swarm onset across sliding time windows of 1 to 12 months. The difference (df = d - d) 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 M6 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. M6 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 Mw ≥ 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 Mw5.9 and Mw5.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 Mw7.3 Hualien earthquake. This swarm-RESs coupling, and its potential connection to the 2024 Mw7.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 deformation 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. 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 Mw7.3 mainshock.

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