Heterogeneous strong asperities and tectonic complexity control irregular cascading ruptures

Relocated aftershocks

We used broadband seismic data from 21 near-field stations (Fig. 1A) to relocate earthquake events between 1 and 28 January 2024, as detailed in Materials and Methods. Using the HypoDD software (21) and conducting four iterations (22), we ultimately identified 13,053 earthquake events. A comparison with the earthquake catalog released by JMA revealed a sharpened spatial clustering in our locations, leading to a more robust depiction of the underlying fault geometries. The aftershock distribution notably differentiates between the western (profile BB′-HH′ in Fig. 1D) and eastern segments (profile II′-LL′ in Fig. 1D) separated at about 137.5°E (Fig. 1C). Most aftershocks in the western section are focused at depths of up to 13 km, following a trend of around 50°. In contrast, aftershocks in the eastern section reach depths of up to 30 km and are spread out along a trend of ~40°. Oppositely dipping trends are observed in profiles HH′ and II′ (~25 km northeast of the epicenter). Profile HH′ and its western counterparts exhibit southeastward dipping bands of seismicity, whereas profile II′ and its eastern counterparts show northwestward dips (Fig. 1). The aftershocks of the 2024 Noto earthquake notably overlap with the 2020–2023 earthquake swarm, indicating that the 2024 Noto earthquake induced rerupture of earthquake swarm–related faults (profile GG′ in Figs. 1D and 2D). We also find that the 2007 M6.9 event likely occurred on different faults from the 2024 M_w 7.5 rupture (profile BB′-CC′ in Fig. 1D), reflecting the complexity of the fault distribution in this area. The substantial overlap between the distribution of 2024 Noto aftershocks and the earthquake swarm, with the hypocenter located in the gap between four clusters of the earthquake swarm, suggests a connection between the nucleation of the 2024 Noto M_w 7.5 earthquake and the earthquake swarms and associated fluids.

Open multimedia modal

A high number of aftershocks is found near 137.25°E, with more than 3000 events around M1 occurring within a 0.05° to 0.1° region (Fig. 2A). Our analysis reveals the presence of a small earthquake swarm triggered by the mainshock, located at a depth of about 13 km (Fig. 2, B and D). This swarm lasted about 1 month, reaching a maximum magnitude of M4.6 (Fig. 2C). The earthquake swarm began on 1 January 2024, with a generally decreasing trend in the frequency of seismic events, effectively concluding around 25 January 2024. Notably, larger seismic events in the swarm area did not generate the usual number of secondary aftershocks, which are commonly seen with typical seismic activity (Fig. 2, B to D). For instance, a M4.6 earthquake on 7 January 2024 occurred during a notable drop in swarm activity.

Slowness-enhanced back projection

Using the SEBP method, we used teleseismic data from the AU and NA arrays to perform a rupture process inversion for the 2024 Noto earthquake (23). Additionally, we used six aftershock events (M > 4.5) to correct for velocity structure heterogeneity, ensuring more robust results (22). The SEBP results effectively illuminate the complex spatiotemporal evolution of the high-frequency radiation of the 2024 Noto earthquake. As shown in Fig. 3, the high-frequency radiations of the M_w 7.5 mainshock indicate that the bilateral rupture mainly occurred in two stages. During the initial 18 s, the rupture propagated slowly southwestward at ~0.6 km/s, showing good consistency in both arrays. Notably, the locations of high-frequency radiations in the first stage overlap well with the earthquake swarm region (Fig. 3A). In the second stage, the rupture began to accelerate bilaterally. From 20 to 40 s, it continued to rupture southwestward at ~2.7 km/s. Between 20 and 56 s, the rupture propagated northeastward at a speed of ~3.1 km/s. Given the resolution limitations of SEBP and the bilateral rupture, particularly in the complex geometry of the swarm-fault region, the apparent back-propagating rupture observed between 0 and 18 s is likely an artifact rather than a real physical process. However, the slow rupture evolution is a reliable feature. It is noteworthy that the strongest high-frequency radiations at the AU array appeared between 30 and 37 s, accompanied by deceleration of the rupture.

Open multimedia modal

Coseismic deformation

We used the pixel offset method on ALOS-2 SAR (Stripmap) data from three tracks to obtain the coseismic line-of-sight deformation field. Additionally, we used differential measurements from 71 GNSS stations within 200 km of the epicenter to obtain three-dimensional coseismic displacements (24). The range pixel-offset maps derived from the ALOS-2 data indicate that the deformation observed in the Noto Peninsula is mainly characterized by uplift, with line-of-sight decreases of up to 4 m at the northwestern end (Fig. 4). The GNSS coseismic displacement field reveals that the nearby stations are predominantly displaced toward the west and uplifting, with a 2-m southwestward displacement near the epicenter and prevailing northwestward movements in the southeastern part of the source region (fig. S1). The general deformation pattern is consistent with a thrust mechanism, coupled with a minor right-lateral component (figs. S1 and S2). The deformation fields obtained from the SAR and GNSS data exhibit a high degree of consistency, indicating the reliability of our observations (Fig. 4, A to C).

Open multimedia modal

Finite fault inversion

Guided by the relocated aftershocks, we performed an inversion using three faults with opposing dips. In the western section, west of 137.5°E, the main model fault dips 45° to the SE. We use a fault that dips 38° to the SE to approximate one or more closely spaced faults that fit the 2020–2023 earthquake swarm. The swarm faults are located beneath the main fault of the 2024 Noto earthquake and are nearly parallel to it (Fig. 2D). While there are apparently multiple subparallel swarm faults, in the inversion, the slip is only inverted for a single planar fault. In the eastern section, the model fault dips 53° to the NW. The hypocenter was placed on the main fault. Our finite fault inversion results fit well with both seismic and geodetic data. The fitting of the horizontal GNSS data is excellent (fig. S1). The SM data are also well fit, but the peak values of simulated ground motions at some stations are slightly lower than the observed values, which may be due to simplifications in the fault and velocity models. Some near-field stations exhibit large errors, possibly because our inversion uses displacement data, and the baseline drift at near-field stations is too large, making it difficult to accurately integrate acceleration into displacement, as seen in station ISK001 (Fig. 4, D to I, and fig. S4). The fitting of the first-order features of the teleseismic waveforms is also good (Fig. 4, J to O, and figs. S5 to S8), although the fitting of high-frequency information is slightly worse, possibly due to limitations in the velocity model (fig. S9). Overall, the simulated and observed data from the optimal model are consistent. A static checkerboard test indicates that the resolution above a depth of ~15 km is high (fig. S10), noting that the slip from the 2024 Noto M_w 7.5 earthquake mainly occurs above 15 km. The geodetic data provide very weak constraints on the eastern fault, so we rely on the SEBP results to analyze the rupture characteristics of the eastern segment. There is certainly some trade-off between slip on the swarm fault and the main fault. However, our resolution test demonstrates that the model and data have sufficient resolution to distinguish between slip on these closely spaced faults and the sequential order of their rupture. This provides a robust foundation for subsequent analyses of the involvement of these faults in the rupture process (fig. S11).

Our preferred finite fault slip model suggests that the total energy released during the 2024 Noto earthquake is equivalent to an event of M_w 7.54. Initially, early bilateral rupture growth occurs within the 2020–2023 earthquake swarm zone, consistent with the findings from the SEBP (Figs. 3 and 5). The initial rupture during the first ~18 s was almost fully confined to within the swarm region. Notably, the initial rupture took place on the main fault of the 2024 Noto earthquake, while the nearby swarm fault started to rupture about 6 s afterward. After 18 s, the slip on the main fault overlapping with the 2020–2023 earthquake swarm area nearly stopped, but the slip on the swarm fault persisted until 36 s. Subsequent rupture of the western segment of the fault involved four main asperities, with the largest one located 40 to 55 km southwest of the epicenter, having a slip value of up to 10 m. The seismic moment (assuming a shear modulus of 30 GPa) of the largest asperity is equivalent to an M_w 7.1 earthquake. The energy release reached its peak between 30 and 35 s, corresponding to the rupture of this largest asperity. Notably, the southwestern rupture also extended to the edge of the asperity located 25 km southwest of the hypocenter at 27 s, leaving a gap of ~10 km by 5 km. Subsequently, this asperity ruptured the gap within 9 s, with the asperity’s shape outlined by the 1-m slip contour (Fig. 5M). During this period, the westernmost asperity ruptured from bottom to top, with a maximum slip exceeding 10 m, accompanied by a surge in the energy release rate of the earthquake (Fig. 5M and figs. S12 and S13). On the eastern segment, there is a shallower asperity above 10-km depth, with a seismic moment equivalent to an M_w 7.0 earthquake. This is consistent with the primary features of the rupture along this segment identified through potency-density tensor inversion results by Okuwaki et al. (7). In the inversion without swarm fault, early slow rupture and delayed asperity rupture were still observed, suggesting the stability of these phenomena (fig. S14). To evaluate the necessity of including swarm fault in the rupture modeling, we compared the two-fault and three-fault models. We found that the three-fault model shows a better SAR data fitting in the swarm fault region (within the black rectangle in Fig. 4C), with the root mean square (RMS) misfit decreased from 14.5 to 12.9 cm, an 11% reduction. This indicates that the swarm fault might have slipped aseismically, eluding detection by seismic data. This provides further evidence for the involvement of earthquake swarm–related fault in the rupture of the 2024 Noto earthquake.

Open multimedia modal

The role of complex swarm faults in initial rupture and aftershock activity

Preexisting complex swarm faults played an important role in the initial slow rupture and post–earthquake aftershock activity. Including the swarm fault in the inversion reduced the SAR data RMS misfit by 11% and relocated aftershocks align along swarm fault (Fig. 2D). By calculating the stress changes on the swarm fault plane induced by either the early slip or the total slip on the main fault, the results consistently indicate that shear, normal, and Coulomb stress changes in the up-dip portion of the swarm fault promote fault slip. This region experiences a Coulomb stress loading of ~1.5 MPa, with a strong spatial correlation between the slip area and the stress-loaded region (fig. S15). This indicates that the swarm fault might have slipped aseismically. These findings support the involvement of swarm fault in the 2024 Noto earthquake rupture. The inversion of the rupture process, using SAR pixel offsets, GNSS, SM, and teleseismic data, reveals that only the main fault experienced rupture in the first 6 s. Following this, the swarm fault began to slip, with both swarm fault and main fault exhibiting bilateral and up-dip rupturing at a very slow speed (~0.6 km/s). This indicates that the rupture of the main fault likely triggered the rupture of the swarm fault, possibly accompanied by aseismic slip. The earthquake swarm activity since 2020 is notably related to fluid influence, revealing that fluid-weakened faults are highly sensitive to stress perturbations and are easily triggered to slip together by nearby fault ruptures. Earlier research also observed slow initial rupture phase (7^, 9^, 10^, 14). In the swarm region, aftershocks of the 2022 M5.1 and 2023 M6.2 events continued migrating upward along the swarm faults after the mainshock, possibly reflecting pore pressure diffusion (4^, 5). Liu et al. (14) provided direct evidence of fluid accumulation through initial rupture analyses of four events in the swarm region. The initial rupture point corresponds to the upper edge of the 2020–2023 earthquake swarm, where fluid accumulation occurred, also implying a possible link between earthquake nucleation and fluid triggering. By using the SEBP technique, we identified strong high-frequency radiators in the earthquake swarm area, which is generally associated with fault strength and geometry changes. In the swarm region, we found that slip on the main fault nearly stopped after 18 s, while slip on the swarm fault persisted for almost 36 s. Resolution tests indicate some uncertainty in the rupture duration of the swarm fault (fig. S11), but the pronounced differences suggest a likely presence of this trend. The slip on the swarm fault likely reflects the combined contributions of one or more small faults within the swarm region, and the extended duration may be attributed to fault weakening. Yoshida et al. (11), on the basis of aftershock relocation, identified several foreshocks occurring on the swarm faults before the mainshock and suggested that rupture on the swarm faults triggered the rupture of the overlying main fault. We propose that smaller foreshocks may have occurred on the swarm fault, which did not exhibit notable slip. The slip on the main fault could have triggered slip on the swarm faults through dynamic and static stress transfer.

Aftershock relocation revealed that the aftershocks from the 2024 M_w 7.5 Noto earthquake were closely aligned with the preexisting faults in the swarm area (Fig. 2D), indicating rerupture of these structures, further underscoring the swarm fault’s role in the earthquake rupture. According to the earthquake catalog provided by JMA, the M_w 7.5 main shock occurred at 16:10:22, but a M5.9 seismic event occurred 2 km west of the main shock 13 s earlier. Due to the overlap of the P-wave arrivals of the main shock and the wave of M5.9 event at near-field stations (fig. S16), we cannot accurately determine the initial time of the M_w 7.5 event. The contribution of the M5.9 event is included in the rupture process inversion and SEBP, which may also be one reason for the unusually slow initial rupture. The 2024 M_w 7.5 earthquake involved multiple faults during its initial rupture phase, with complex fault geometry, fluid-induced fault weakening, and interactions among multiple faults being crucial factors in the slow rupture observed in the first 18 s.

Delayed rupture behavior of strong asperity

Our finite fault inversion findings indicate a fascinating rupture process involving delayed failure of an asperity during the 2024 Noto earthquake. Before 27 s, the southwestward rupture left a gap measuring 10 km by 5 km, located about 25 km southwest of the epicenter. Between 27 and 36 s, the asperity within this gap experienced a delayed rupture. Initially, the robust asperity notably impeded the slip, causing the rupture to navigate around it from 18 to 27 s after the rupture began. During this time, the strong asperity acted as a barrier, but the slip occurring around it increased stress on the asperity. Eventually, the dynamic and static stress loading overcame the rock strength, resulting in rapid failure of the strong asperity. The delayed rupture of strong asperities can be a secondary source of powerful ground motions, sometimes referred to as compound rupture (25). Xu et al. (10) also found another delayed rupture asperity located 10 km southwest of the epicenter, to the east of the asperity identified in this study. We believe that there may be more than one instance of compound rupture behavior in that heterogeneous area, but the different fault geometries used have varying levels of analytical capability. Numerical simulations have indicated that increasing local fault strength can lead to such delayed ruptures (26). This compound model suggests that barriers not only can arrest ruptures but also can create additional sources, which can exacerbate earthquake damage.

We found that the rupture mechanism of the peak-slip asperity at the far western end was a bottom-to-top rupture instead of a lateral rupture in the horizontal direction of propagation (fig. S13). This could again be due to stronger locking at the shallower asperity, with the shallow asperity following along after slip had already occurred below. This can also be referred to as a compound rupture.

Cascading rupture mechanism in complex earthquake swarm region

Elucidating the mechanisms of large cascading ruptures in earthquake swarm regions allows for a better understanding of the large-scale seismic failure of complex fracture networks. Earthquake swarms are often linked to the movement of fluids and/or aseismic slip and are generally characterized by small to moderate earthquakes, with only a few experiencing larger events. For instance, the largest earthquake of the Matsushiro earthquake swarm in central Japan, lasting over 5 years, was a M5.4 event, and the L’Aquila earthquake swarm in Italy had a seismic event of M_w 6.3 (27^, 28). However, earthquakes larger than M7 are rare in earthquake swarm areas. To enhance our understanding of the mechanisms behind this large earthquake, we conducted a comprehensive analysis of fault structures and rupture processes of the 2024 Noto earthquake event. First, both finite fault inversion and SEBP results indicate that the earthquake rupture is mainly divided into two stages. The first phase is largely confined to the 2020–2023 earthquake swarm area, during which the rupture on the main fault triggered the rupture of underlying swarm fault segments. This prolonged slow rupture includes the contribution from the M5.9 foreshock and may also involve the M5.9 foreshock triggering the onset of the M_w 7.5 main shock through dynamic stress. After 18 s, the event transitioned into the second stage. The slow rupture further triggers the main asperity in the southwest direction by stress changes, also triggering another fault in the northeast direction, representing an M_w 7.0 subevent. This forms a typical cascading rupture event involving multiple faults, which is observed in several large earthquakes, such as the 2016 M_w 7.8 Kaikoura Earthquake and the 2021 M_w 7.4 Maduo earthquake (29^, 30). In summary, complex fault structures and inhomogeneous asperities played a crucial role in the extensive cascading rupture of the M_w 7.5 Noto earthquake in 2024.

Characteristics of the 2024 earthquake swarm and its distinct aftershock patterns

We find that the 2024 M_w 7.5 Noto mainshock triggered a small earthquake swarm located 3 km above the 2020–2023 swarm (Fig. 2, B and D). We assumed the receiver fault geometry of the modeled swarm fault and adopted a generic effective coefficient of friction of 0.4 to calculate the stress changes on the swarm fault. Although the normal stress change derived from the mainshock’s slip distribution is compressive, the shear stress loading is relatively high. As a result, the Coulomb stress change indicates that stress loading exceeding 1 MPa occurred at the location of the 2024 swarm (see Materials and Methods and fig. S17), providing direct evidence that the 2024 Noto mainshock statically triggered the small earthquake swarm. This recent swarm began on 1 January 2024 and subsided around 25 January. After 6 January, the number of earthquakes started to decay. Notably, the larger quakes in this swarm did not produce many aftershocks, which is unusual compared to typical earthquake sequences. Until 25 January, the swarm events released energy comparable to an M_w 5.4 earthquake. Figure 2C shows the earthquake frequency in an area of 0.1° × 0.1° within the 2024 swarm region. Two notable events (M4.6 on 7 January and M4.5 on 12 January) occurred during declines in earthquake rate, and the event rate continued to drop afterward. To verify whether larger events do not produce a rich aftershock sequence due to noise levels, station outages, or during times of frequent events, we used the maximum curvature method to calculate the daily changes in the minimum completeness magnitude (M_c) for the 28 days following the M_w 7.5 Noto earthquake in 2024 (31). Our findings indicate that, aside from the first 3 days—when a higher number of seismic events during early aftershock activity resulted in a larger M_c—M_c remained relatively stable around M1 for the remainder of the period (figs. S18 and S19). This suggests that the phenomenon of larger seismic events occurring without a rich aftershock sequence is a valid observation. A similar pattern was noted in the 2020–2023 swarm (Fig. 2E), where many events greater than M3.5 also did not lead to numerous aftershocks. Lei et al. (32) used the epidemic-type aftershock sequence (ETAS) model and found that earthquakes caused by hydraulic fracturing do not produce as many aftershocks as tectonic earthquakes do. This may be due to the absence of preexisting fault damage zones in fluid-induced earthquakes, which are typically found in tectonic events. Consequently, fluid-induced seismic events may not trigger aftershocks in the same way as conventional earthquakes.

Aftershock relocation

To compile a precise catalog of aftershocks following the 2024 Noto earthquake, seismic data were gathered from 21 broadband stations with an epicentral distance of 100 km. These data were obtained from Hi-net, F-net, and JMA networks operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). The data collection period spanned from 1 to 28 January 2024. The aftershock zone was divided into a western segment beneath the Noto Peninsula with dense station coverage and an eastern segment offshore with less station coverage, primarily concentrated on the landward side in the SE direction. Notably, a station was also positioned on an island to the NW. The seismic waveform data were processed using the PALM code through a two-step procedure (22). First, a joint algorithm using STA/LTA and kurtosis was used to identify pairs of P and S arrivals. This involved band-pass filtering continuous waveforms between 1 and 20 Hz, followed by a search for S waves within a 16-s window after the P arrivals. Second, a matched-filter detection technique was applied with a cross-correlation threshold of 0.2. For differential travel time calculation, we used a cross-correlation threshold of 0.3. To ensure accurate association of the same seismic event, a minimum of four stations with an original time deviation of less than 2.0 s and a P travel time residual of less than 1.5 s were required. We adopted the one-dimensional (1D) velocity model by JMA for the analysis (fig. S9).

Slowness-enhanced back projection

Back projection is a method of array processing that relies on the similarity of seismic waveforms within the same array. It can elucidate the kinematic characteristics of rupture, such as rupture direction, velocity, and energy release size, with minimal prior knowledge of fault geometry or rupture velocity. Traditional back projection assumes a uniform velocity structure within the source region, which may lead to errors in cases of large and laterally heterogeneous source regions. To address this limitation, the SEBP method is used, which incorporates spatially varying slowness correction terms derived from aftershock events to better account for the heterogeneous structure and provide more accurate source process information (23).

The AU array to the southwest and the NA array to the northeast of the study area were selected because of their good constraints on the rupture process. Teleseismic distances between 30° and 90° were chosen to mitigate the influence of strong depth phases. A total of 121 and 322 broadband seismograms from the AU and NA arrays, obtained from the Incorporated Research Institutions for Seismology (IRIS), were separately analyzed (fig. S20). The Multiple Signal Classification back projection method was applied, aligning P arrivals through band-pass filtering across multiple frequency bands (33). Six aftershocks with magnitudes above M4.5 distributed throughout the aftershock zone were selected to evaluate the accuracy of the velocity structure.

We found that the average spatial bias of back projection–inferred locations for the AU and NA arrays, based on a uniform velocity structure, was 11.1 and 12.1 km, respectively. To mitigate spatial bias caused by crustal heterogeneity, slowness correction terms were derived from the spatial bias of the six selected aftershock events (tables S1 and S2). After applying spatially variable slowness corrections across the entire source region, the average spatial bias improved notably to 5.2 km for the AU array and 7.7 km for the NA array (fig. S21). This enhanced accuracy allows for a more detailed understanding of the spatiotemporal evolution process of the earthquake rupture.

SAR pixel offsets and GNSS observations

We used SAR data obtained from the L-band ALOS-2 satellite operated by the Japan Aerospace Exploration Agency. Three pairs of SAR images were processed from ascending tracks 121 and 127, as well as descending track 26, as outlined in table S3. We processed the data to obtain pixel range offset maps of the 2024 Noto earthquake through the pixel offset-tracking (POT) method (34). Due to the considerable coseismic deformation causing severe decorrelation, a matching window size of 120 pixels in range and 84 pixels in azimuth was used, with a moving step of 30 × 21 pixels. The POT method was implemented on each of the three consecutive frames, and the resulting frames were combined to create a comprehensive deformation map. To improve efficiency in the finite fault inversion, a quadtree method was used to downsample the deformation field, resulting in 944, 701, and 721 points for inversion (35).

Japan maintains a dense GNSS network operated by Geospatial Information Authority of Japan. To study coseismic deformation, 71 stations located within 200 km of the earthquake epicenter were chosen on the basis of data quality and proximity. Data from these stations were collected at 30-s intervals. GNSS data recorded on 1 January 2024, underwent processing using precise point positioning with ambiguity resolution software (36). This processing involved using final orbits from the Center for Orbit Determination in Europe, phase clocks reestimated by Wuhan University, and phase bias products, with a cutoff elevation angle of 10°. Thirty-second time series on a single day were generated for the 71 selected stations in the International Terrestrial Reference Frame 2014. Following established methodologies, coseismic displacements were determined by subtracting the average preseismic position from the postseismic position within a 50-min time frame before and after the rupture (fig. S22) (37^, 38).

Local SM and far-field teleseismic data

Japan maintains a closely spaced network of SM stations, known as k-net and kik-net, overseen by NIED. We identified 18 SM stations located within a 150-km radius. These stations, unlike broadband seismometers, do not experience instrument clipping and offer high-frequency sampling rates, enabling detailed examination of rupture processes (25^, 39). For inversion efficiency, acceleration data from SM stations typically need to be integrated twice to derive displacement time series for rupture process inversion. However, SM data often suffer from substantial baseline drift. To mitigate this effect, we applied band-pass filtering in the frequency range of 0.02 to 0.125 Hz. Nevertheless, for some near-field stations with severe baseline drift, integration to displacement may still introduce certain errors, as observed at station ISK001 (fig. S4).

We chose broadband teleseismic waveforms from 48 stations at epicentral distances spanning 30° to 90° sourced from the IRIS database, ensuring a comprehensive azimuthal distribution. Subsequent to the removal of instrument response, the body waves (P and SH waves) were subjected to bandpass filtering with corner frequencies set at 0.005 to 1 Hz, while the surface waves (Rayleigh and Love waves) underwent similar bandpass filtering with corner frequencies ranging from 0.005 to 1 Hz.

Inversion method and model parameters

We used the methodology introduced by Ji et al. (40) to conduct a joint inversion of SAR pixel offsets, static GNSS, SM, and teleseismic data for the detailed rupture characterization of the 2024 Noto M_w 7.5 earthquake. This approach involves nonlinear inversion techniques to optimize the agreement between observed and simulated data to derive parameters such as slip amplitude, rake, rupture time, and duration for individual subfaults (41). The optimization process is carried out using simulated annealing. A source time function was constructed for each subfault using an asymmetric cosine function with eight half-durations of 1.5 s to account for variable rupture velocity and duration. The configuration of fault geometry was determined from the relocated aftershocks. The western segment was defined as centered on the epicenter, extending 98 km along a 58° strike and 35 km along a 45° dip angle, with the upper edge of the fault positioned close to the surface. The eastern segment extends 55 km along a 40° strike and 30 km along a 53° dip angle. To enhance our understanding of the relationship between the initial rupture, the main fault, and the faults associated with the earthquake swarm, we incorporated the swarm fault into the inversion analysis. By analyzing the earthquake swarm from 2020 to 2023, we found that two main fault planes are situated closely together in a vertical alignment, which are both positioned beneath the main fault (Fig. 2D). Accordingly, our inversion uses a single fault plane to represent the contribution of one or more swarm-related faults (Fig. 3A). The fault was discretized into subfaults of dimensions 3.5 km by 3.4 km. For the inversion of the regional data, we adopted the 1D velocity model proposed by JMA and used the frequency–wave number integration method to compute Green’s functions (42). The propagator matrix approach was used for teleseismic data to calculate Green’s functions (43). Laplacian regularization was applied over the fault plane to ensure the spatiotemporal stability of slip evolution. We set the rupture velocity range to 0.3 to 3.5 km/s, with a reference velocity of 2.5 km/s. We use moment magnitude as the target in the inversion process, automatically searching for an appropriate smoothing parameter that matches the target magnitude. Different observational data have different uncertainties. Given that SH waves have a lower signal-to-noise ratio and greater arrival time errors, we assign P waves a weight twice that of SH waves. The Green’s function of long-period surface waves is more reliable than that of body waves (44), so we assign surface waves a weight twice that of P waves. Because near-field data provide higher resolution, we give SAR pixel offsets, horizontal GNSS, and SM data the same weight as surface waves. Vertical GNSS, which is less accurate than horizontal GNSS, is assigned one-quarter the weight of horizontal components (45). Thus, the data are weighted as follows: SAR pixel offsets:horizontal GNSS:vertical GNSS:SM:P:SH:Rayleigh:Love = 4:4:1:4:2:1:4:4.

Resolution tests

To assess the resolution capability of the data and model in capturing rupture evolution, we conducted three resolution tests using the aforementioned parameters. First, static data (SAR, GNSS, and static SM offsets) were used to evaluate the ability of the data and model to constrain the static slip distribution (fig. S10). Second, to investigate the resolution capability of the dataset in distinguishing slip on two closely spaced parallel faults (main fault and swarm fault), we performed a resolution test using all available datasets (SAR, GNSS, SM, and teleseismic data). In the synthetic test, two 6.5-m slip asperities were assigned to the main and swarm fault, respectively. The asperity on the main fault ruptured between 5 and 11 s, while the asperity on the swarm fault ruptured between 13 and 17 s. In the inversion, we imposed two separate nucleation points on the main fault and the swarm fault. This test aimed to evaluate the ability of the data and model to resolve the sequential rupture of the two faults (fig. S11). Third, to evaluate the necessity of including a swarm fault in the rupture inversion, we performed an additional rupture inversion excluding the swarm fault while using the same dataset as the three-fault inversion (SAR, GNSS, SM, and teleseismic data) (fig. S14).

Stress change calculation

We calculated the shear stress, normal stress, and Coulomb stress changes on the swarm fault using the slip distribution on the main fault to evaluate the triggering effect of the rupture on the swarm fault. The stress changes of each subfault were computed at its center to evaluate the distribution of stress changes encouraging slip in the coseismic slip direction. Snapshots of the rupture process indicate that notable slip on the swarm fault occurred after 18 s (Fig. 5). Therefore, we estimated the stress changes imparted to the swarm fault by the slip on the main fault within the first 18 s. To further assess whether the stress changes remained encouraging and comparably high after the rupture had completed, we also calculated the stress changes on the swarm fault induced by the total slip on the main fault. The rake angle of each subfault was set according to its coseismic slip vector (strike, 51; dip, 38; and rake, 120). We assume an effective coefficient of friction of 0.4. The results consistently reveal a relatively high Coulomb stress loading in the up-dip direction of the swarm fault (fig. S15).

To assess the triggering effect of the 2024 Noto mainshock on the 2024 earthquake swarm, we calculated shear stress, normal stress, and Coulomb stress changes induced by the mainshock slip. As shown in Fig. 2D, the triggered swarm events are primarily concentrated at depth of 13 km. Therefore, we computed stress changes at depth of 13 km using the slip distributions of the three faults. Given the proximity of the 2024 earthquake swarm to the swarm fault and the possibility of rupture on a smaller fault within the swarm fault, the receiver fault was set with the swarm fault parameters (strike, 51°; dip, 38°; rake, 120°) and adopted a generic effective coefficient of friction of 0.4. The results indicate that the 2024 earthquake swarm at depths of 13 km is located within the Coulomb stress loading region (fig. S17).

Acknowledgments