A Multiplex Rupture Sequence Under Complex Fault Network Due To Preceding Earthquake Swarms During the 2024 Mw 7.5 Noto Peninsula, Japan, Earthquake

A devastating earthquake with moment magnitude 7.5 occurred in the Noto Peninsula in central Japan on 1 January 2024. We estimate the rupture evolution of this earthquake from teleseismic P‐wave data using the potency‐density tensor inversion method, which provides information on the spatiotemporal slip distribution including fault orientations. The results show a long and quiet initial rupture phase that overlaps with regions of preceding earthquake swarms and associated aseismic deformation. The following three major rupture episodes evolve on segmented, differently oriented faults bounded by the initial rupture region. The irregular initial rupture process followed by the multi‐scale rupture growth is considered to be controlled by the preceding seismic and aseismic processes and the geometric complexity of the fault system. Such a discrete rupture scenario, including the triggering of an isolated fault rupture, adds critical inputs on the assessment of strong ground motion and associated damages for future earthquakes.


Introduction
Growing observational evidence has unveiled diverse earthquake rupture behaviors, which showcase for example, multiplex triggering ruptures (Meng et al., 2012;Hicks & Rietbrock, 2015;Fan et al., 2016;Vasyura-Bathke et al., 2024), transient rupture deceleration and acceleration (Bao et al., 2019;Ding et al., 2023), or apparent backward rupture propagation (Hicks et al., 2020;Vallée et al., 2023).The configuration of a complex fault network can be a key driver of multi-scale cascading ruptures with alternate rupture directions (Yamashita, Yagi, & Okuwaki, 2022;Ohara et al., 2023;Okuwaki et al., 2023).Even in fault systems with apparently simple geometries, complicated bypassing or boomerang-like rupturing can be induced by inhomogeneous barriers on fault surfaces due to heterogeneous fracture energy or strength (Hicks et al., 2020;Yagi et al., 2024).Complex rupture sequences on differently oriented faults can impact the generation of strong ground motions (e.g.Aochi & Madariaga, 2003;Chu et al., 2021;Yagi et al., 2023;Taufiqurrahman et al., 2023).Robust estimates of diverse rupture behaviors benefit our understanding of earthquake-source physics and provide critical inputs for better assessments of earthquake hazards, especially those involving a complex fault network within heterogeneous media.
The Noto Peninsula in central Japan (Fig. 1) is one such environment.Active faults are distributed in the off-shore region along the northern coast of the peninsula, forming a segmented fault network (AIST, 2012).The active faults are characterized by reverse faulting that has resulted from northwest-southeast-oriented compression; these faults are considered to have been re-activated by inversion tectonics after their formation under normal faulting during the opening of the Japan Sea before 15-20 Ma (e.g.Sato, 1994;Okamura et al., 1995).It is not rare to have modest to large magnitude (M) 6 class earthquakes in this region; recent major events include the 2007 M 6.9 and the 2023 M 6.5 earthquakes (Japan Meteorological Agency, 2024;Hiramatsu et al., 2008;Yoshida, Uno, et al., 2023) (Fig. 1).The 2023 M 6.5 event, which occurred in the northern tip of the peninsula, was preceded by the long lived earthquake swarms from early November 2020 (Amezawa et al., 2023;Yoshida, Uno, et al., 2023;Yoshida, Uchida, et al., 2023;Kato, 2024;Shelly, 2024).Geodetic analyses suggest a series of transient aseismic to seismic deformation from deep to shallow domains, that is aseismic deformation caused by fluid spreading from depth, which triggered up-dip earthquake swarm activity (Nishimura et al., 2023).Thus, the Noto Peninsula is an intriguing environment in which earthquake ruptures can occur on geometrically complex fault segments and in spatially heterogeneous strength fields due to the existence of fluids.
On 1 January 2024, a moment magnitude (M W ) 7.5 earthquake occurred in northern Noto Peninsula (Japan Meteorological Agency, 2024) (Fig. 1).The 2024 Noto Peninsula earthquake severely affected people and infrastructure due to strong shaking, causing at least 245 causalities as of 16 May 2024 (Fire and Disaster Management Agency, 2024).The global centroid moment tensor project (Dziewonski et al., 1981;Ekström et al., 2012) finds the earthquake is characterized by reverse faulting due to northwest-southeastoriented compression.Aftershocks determined by the Japan Meteorological Agency (Japan Meteorological Agency, 2024) extend out about 150 km in the northeast-southwest direction, slightly changing strike with a subtle S-shaped curve (Fig. 1).The dip directions of the aftershocks are not obvious, but they roughly show southeastward dips on the western side of the epicenter and northwestward dips on the eastern side (Fig. S1).Synthetic Aperture Radar (SAR) data (Geospatial Information Authority of Japan, 2024) show two major patches of deformation on the western side of the epicenter, including vertical displacement of up to about 4 m.A tsunami was also recorded at gauges along the coast (Fujii & Satake, 2024;Mizutani et al., 2024;Kutschera et al., 2024).
These early observations collectively suggest that the mainshock fault geometry is not as simple as one defined by a solely rectangular plane; rather, it appears to include curved and/or segmented planes.It is inherently challenging to understand the causal relationship between rupture evolution and the geometric features of faults on an observational basis, because inappropriate assumptions about fault geometry can be a major source of modeling error and easily bias the solution, making robust interpretation of source processes difficult (e.g.Ragon et al., 2018;Shimizu et al., 2020).For our analysis of the 2024 Noto Peninsula earthquake, we apply a recently developed potency-density tensor inversion (PDTI) method that can flexibly estimate rupture evolution without strong assumptions about fault geometry (Yagi & Fukahata, 2011;Shimizu et al., 2020;Yamashita, Yagi, Okuwaki, Shimizu, et al., 2022).We discuss a multiplex rupture sequence across the segmented fault network and its relation to the preceding earthquake swarm, which together control the nucleation process and irregular rupture behavior of the 2024 Noto Peninsula earthquake.

Data and Methods
Inappropriate assumptions about fault geometry are a major source of modeling error in finite-fault inversions (Ragon et al., 2018;Shimizu et al., 2020;Dutta et al., 2021).A realistic fault geometry can be assumed by prescribing one or multiple finite-fault plane(s) guided by active faults, slab geometry models, or aftershock distributions.However, this is not necessarily the true fault geometry hosting the co-seismic rupture of interest.In particular, the conventional finite-fault inversion method requires that the fault slip vectors be strictly forced to span the prescribed model surface, which is a strong constraint that may induce modeling errors when the fault geometry assumption is inappropriate, leading to biases in both the modeling and interpretation.A recently developed PDTI method (Shimizu et al., 2020) solves for the potency-rate density tensor distribution by alternatively representing a fault slip as a composite of five-basis moment tensors (e.g.Kikuchi & Kanamori, 1991).This composite is independent of the prescribed model fault plane and thus enables the simultaneous estimation of both the rupture evolution and fault geometry.This flexible approach helps to avoid modeling bias and leads to robust modeling of earthquake source processes, and it has been proven to be efficient in modeling diverse rupture processes for a variety of tectonic environments or unknown fault configurations (e.g.Shimizu et al., 2020;Tadapansawut et al., 2021;Okuwaki & Fan, 2022;Yagi et al., 2023;Ohara et al., 2024).
For the 2024 Noto Peninsula earthquake, both the active faults and the aftershocks show spatial changes in the fault strike and dip orientations (Figs. 1 and S1).Such features motivate us to use the PDTI method to analyze the source process of the earthquake.It is important to reiterate that this PDTI approach does not require a strict assumption about the fault geometry, and the adopted model space geometry (model fault) does not necessarily coincide with the actual fault geometry.Such a flexible representation of earthquake sources can be made by projection of an intrinsically volumetric distribution of potency-rate density tensors onto a 2-D model fault (Shimizu et al., 2020).Teleseismic body waves have rich information about source evolution, including the variation of focal mechanisms (Kikuchi & Kanamori, 1991;Shimizu et al., 2020Shimizu et al., , 2021)), although the spatial resolution is relatively low compared to other near-field datasets (e.g.SAR).The PDTI method using teleseismic body waves is not very sensitive to source location errors, which is in turn a theoretical ground that our PDTI approach does not require a strict fault geometry assumption.However, the dip angle of the model geometry may affect the final solution, because the Green's functions change with depth.Therefore, the horizontal location of slip can vary according to the horizontal location and dip direction of the model faults.We build source models based on two alternative model faults (Fig. S1): (1) a twin model fault with a composite of southeast-and northwest-dipping planes (strike/dip: 55 • /35 • and 225 • /35 • , respectively) and (2) a single model fault with a southeast-dipping plane (55 • /35 • ).For the twin model fault, fault lengths are 95 km in the western section and 55 km in the eastern section.For the single model fault, the fault length is 150 km.For both models, the fault width is 35 km, and the model faults are discretized into 5 km × 5 km subfaults along the strike and dip directions.In the results and discussion, we primarily focus on the results of the twin model fault, but we also discuss how this model setup may affect the solution and our interpretation of the source process of the 2024 Noto Peninsula earthquake.
We adopt a hypothetical model rupture front that defines a timing of rupture initiation at each source element, propagating at 3.9 km/s based on the S-wave velocity around the source region (Table S1), which is fast enough to allow for various rupture scenarios.We represent a slip time function at each source element by a series of linear B splines at 0.6-s intervals.We set the duration of the slip time function at each source element from the hypothetical rupture initiation to the rupture termination, which is set 45 s after the start of the earthquake.Such a wide model space enables us to flexibly represent potentially irregular rupture or slip behaviors, including multiple peak of slips, reverse migration of the rupture front, or re-rupture at a given source point.
We use the vertical component of the teleseismic waveforms at 32 globally operated broadband stations via the SAGE Data Management Center (Fig. 1c).The data are selected based on clear P -wave arrivals that can be reliably picked and to achieve good azimuthal coverage (e.g.Okuwaki et al., 2016).The data are then resampled at 0.6 s intervals after removing instrumental responses and converted into velocities.Green's functions are calculated based on the method of Kikuchi and Kanamori (1991).We adopt the ak135 velocity structure model (Kennett et al., 1995) for the near-source region to calculate the Green's functions.Our selection of the velocity model is evaluated by also adopting the CRUST 1.0 model (Laske et al., 2013) (Table S2), although we find that the selection of the velocity model does not significantly affect the solution (Fig. S2).The JMA registered the mainshock as having a JMA magnitude (M JMA ) 7.6 occurring at 37.496 • N, 137.271 • E at 2024-01-01 07:10:22 (UTC).However, 13 s earlier (at 2024-01-01 07:10:09 UTC), they also register a M JMA 5.9 foreshock at 37.508 • N, 137.230 • E (Japan Meteorological Agency, 2024).After inspecting the teleseismic dataset, clearer P -wave onsets can rather be identified for the signals associated with the preceding event.Thus, we regard both the M JMA 5.9 and 7.6 events as a series of the mainshock sequence, and we adopt the initial rupture point (hypocenter) to be at 37.508 • N, 137.230 • E, and 12 km depth for our modeling.

Results
As shown in Figs.2-4, the source model built by using the PDTI method for the 2024 Noto Peninsula earthquake shows a series of discrete rupture episodes (E1, E2, and E3) following a long (about 10 s) and quiet initial rupture episode (E0) around the hypocenter.The timing and location of these rupture episodes are shown in Fig. 4. The total seismic moment is 2.5×10 20 N m (M W 7.5).The general faulting mechanism extracted from the space-time integration of all the potency-rate density tensors is characterized by reverse faulting with northwest-southeast compression (Fig. 1), but it involves considerable changes of fault geometry in space and time, which are described below.In this section, we primarily focus on the results built upon the twin model fault, although both models of the single and twin model faults share common features of rupture evolution and fault geometry with comparable data fits (Figs. S7 and S8).Differences between the models that may affect the interpretation of the results are discussed in Section 4.

A quiet initial rupture episode: E0
During the first 10 s after the origin time, there is a minor initial rupture episode around the hypocenter (E0), which corresponds to 1% of the total seismic moment (M W 6.3).The best-fitting double couple solution extracted from the resultant potency-rate density of this episode shows reverse faulting with a strike of 68 • at the maximum location of the potency rate.The rupture direction is not clearly inferred, but it mainly propagates west from the hypocenter (Fig. 2).

A minor rupture episode: E1
Ten seconds after the origin time, a higher potency-rate event than E0 begins around from 0 km to 40 km west of the hypocenter; this episode (E0) continues until 24 s.The strike angle extracted from the potency-rate density tensors of this episode is 55 • at the maximum location of the potency rate, which is the same as the strike angle of the corresponding section of the model fault (Figs. 3 and 4).The rupture of this event is mostly shallow (< 15 km).
3.3 A major rupture episode in the western section: E2 In the western section of the model fault, a major event (E2) occurs from 25 to 36 s after the origin time.It extends from 0-65 km west of the hypocenter.The strike angles during the episode rotate counter-clockwise from those in E1, toward the north-south direction, with a range of 21 • to 48 • (Fig. 4).The strike at the maximum location of the potency rate is 37 • .This event ruptures down to 20 km depth, which is deeper than E1.

A major rupture episode in the eastern section: E3
At almost the same time as E2, another major rupture episode (E3) occurs in the eastern section during 26-34 s after the origin time.It extends 10-60 km east of the hypocenter.The strike angles during E3 also rotate counter-clockwise relative to those in E1, with a range of 215 • to 224 • , which are more aligned than during E2 (Fig. 4); the strike at the maximum location of the potency rate is 215 • .Note that for E3, the strike angle is selected from the ones of the double-couple solution that is close to the northwest-dipping model plane.As discussed in the later section, the highest potency rate is generally obtained in the shallow part for both the twin and single fault models; because of the alternate dipping sense of the model faults, the twin fault model finds the highest potency rate close to the southern side of the model fault, whilst the single model finds it to the northern side of the model fault (Figs. 4,S4 and S6).The rupture suddenly becomes faint after 34 s, with a minor moment-rate release from 38 to 42 s.

Discrete rupture episodes robustly estimated against model setups
The source process of the 2024 Noto Peninsula earthquake can be characterized by multiple spatiotemporally segmented rupture episodes.The rupture episodes in the western (E1, E2) and eastern (E3) model domains are separated by the hypocenter region, which only hosts the initial minor rupture episode (E0).To evaluate whether this discrete rupture process is artificially made by the segmented domains designed for the twin model fault, we perform the same inversion procedure, but using an alternative model-fault geometry composed of a single southeast-dipping plane without imposing any prescribed segments (Figs.S3-S5).Even in this alternative case, the subtle initial rupture episode (E0) is still robustly estimated, followed by the E1 and E2 rupture episodes in the western side.Near-field records from local stations show a relatively larger amplitude at a station on the western side of the epicenter (Fig. S11), which may indicate the rupture directivity during the first 10 s and is consistent with the westward rupture direction of E0.E3 is also found in the eastern section across the hypocenter.The strike changes among the rupture episodes are similar to those estimated for the twin fault model; the strike orientations of both E2 and E3 rotate counter-clockwise from the reference model strike (55 • ) (Fig. 4).These common features retrieved from the two different model setups suggest that the segmented rupture episodes separated by the hypocenter are robust features of the 2024 Noto Peninsula earthquake.
A notable difference of the results between the two models is their slip locations in the eastern section.For the twin fault model adopting the northwest-dipping segment in the east, the highest potency rate is located in the south-eastern side of the eastern fault, whereas it is in the north-western side in the single fault model (Fig. S6).In this case, our modeling approach leaves an ambiguity of the dipping sense, mostly due to the lower spatial resolution of the teleseismic data, as compared to that of the near-field data sets.Absolute slip locations will need to be further evaluated with other data sets that have enough sensitivity to slip location (e.g. with tsunami data), but this is beyond the scope of the present study.That said, our approach is advantageous to robustly estimate rupture evolution and fault geometry changes against a selection of model setup.In addition, this exercise adopting conjugate dipping planes suggest that our modeling approach can constrain rupture depth, which is concentrated in the shallower domain for both model faults.

Control of rupture behavior due to preceding earthquake swarms
Our source model shows the initial rupture episode (E0) is quiet and lasts 10 s before the main rupture begins.As mentioned previously, we regard the preceding M JMA 5.9 and the subsequent M JMA 7.6 events as one continuous earthquake.The near-field strong motion records show that the M JMA 5.9 event had an unusually long duration, when compared with other similarly sized events, and did not terminate before the M JMA 7.6 event initiated (Fig. S11).Irregular initial rupture phases have been observed for other earthquakes (e.g.Ellsworth & Beroza, 1995;Abercrombie & Mori, 1994), although the 10-s duration of our E0 episode is longer than most cases documented in the literature.One exception is the 2014 M W 8.1 Iquique, Chile, earthquake, which had a long (20 s) initial rupture phase (e.g.Yagi et al., 2014) and was preceded by swarm-like foreshock activity that lasted around 2 weeks (e.g.Ruiz et al., 2014;Kato & Nakagawa, 2014).These observations raise the question of what environment fosters the development of the type of rupture behavior observed in the 2024 Noto Peninsula earthquake.
The northern tip of the Noto Peninsula near the hypocenter of the 2024 Noto Peninsula earthquake has experienced long-lived earthquake swarm activity since early November 2020, which seemingly ended when the adjacent M 6.5 earthquake occurred on 5 May 2023 (Amezawa et al., 2023;Kato, 2024) (Figs. 1 and S9).Geodetic studies have found that aseismic crustal deformation accompanies the swarm activity, suggesting upwelling fluid migration from around 16 km depth to the shallow permeable fault zone, which triggers a series of diffusive aseismic processes to seismic deformation (Nishimura et al., 2023).The E0 domain overlaps the earthquake swarm region, and the eastern edge of E1 coincides with the western edge of the swarm region (Figs. 3 and 4).It is likely that the E0 domain experienced strain release preceding the mainshock due to the fluid-induced earthquake swarms and associated aseismic deformation, making the local environment unfavorable to a spontaneous fast rupture upon a large stress drop.This development prior to the mainshock is considered to be the cause of the long-lasting initial rupture phase of E0.Earthquake-swarm activity has often been interpreted to be associated with aseismic slip and/or fluid migration (e.g.Vidale & Shearer, 2006;Shelly et al., 2016;Fukuda, 2018;Ross et al., 2020;Nishikawa et al., 2021;Im & Avouac, 2023), and has been evaluated in numerical simulations under variable physical conditions (e.g.Zhu et al., 2020;Dublanchet & De Barros, 2021;Wang & Barbot, 2023).For example, a model incorporating permeability changes has predicted fluid-driven aseismic slip and swarm seismicity, together with fluid pressurization ascending through the seismogenic zone, which can eventually lead to a large earthquake (Zhu et al., 2020).Here the 2024 Noto Peninsula earthquake provides additional observational evidence that the co-seismic rupture process itself is also controlled by the preceding fluid-induced swarms, which contributes to the irregular initial rupture phase and the following multiplex rupture episodes across the swarm region (Figs.2-4).

Triggering of segmented fault ruptures
The strike orientations differ among the rupture episodes.They are relatively eastwest for E1 in the central section and relatively north-south for E2 and E3 in the western and eastern sections.Such differences are robustly retrieved when adopting the alternative model fault geometry (Fig. 4).In addition, the changes of strike are also consistent with the aftershock distribution (Japan Meteorological Agency, 2024) and the surface displacement pattern observed inland by SAR (Geospatial Information Authority of Japan, 2024).The distinct features of these strike orientations suggest that E1 and E2 ruptures differently oriented fault segments, evolving from the smaller-scale rupture of E1 to the larger-scale rupture of E2.Such multi-scale rupture growth controlled by geometric features of the fault system is similarly observed during the 2023 Türkiye and Syria earthquake doublet (e.g.Okuwaki et al., 2023).The duration of the western rupture (E1 + E2) is longer than that of the eastern rupture (E3) (Fig. 2b), and the strike orientations in the western section display a wider range of variability (Figs. 4 and S10).Although for convenience, we have described our results by defining rupture episodes (e.g.E2), given the considerable variability of the fault geometry seen during E2, the E2 domain is not necessarily a single fault.It is more likely to comprise multiple faults that take a longer time to be ruptured, and may be related to the complex aftershock distribution in the western side of the peninsula (Fig. 1).
The E3 rupture is isolated from the other episodes, and a rupturing path either from the hypocenter or the other episodes is not obvious from our solution.Although it is generally difficult to rigorously estimate rupture speed from our solution due to the limited spatial resolution and smoothing effects, the apparent migration speed from E1 to E3 is roughly estimated to be ∼3.4 km/s, which is close to the S-wave velocity around the source region (e.g.Table S1).Distant rupture episodes such as E1 possibly contribute to dynamic disturbances in the E3 domain and the resultant rupture.As discussed in the previous subsection, the absolute location of slip and the dip direction of E3 are difficult to uniquely determine solely by our modeling using teleseismic data, but the northwestdipping E3 rupture found with the twin fault model likely corresponds to one of the potential tsunami source faults recognized by the Ministry of Land, Infrastructure, Transport and Tourism and Japan Sea Earthquake and Tsunami Research Project before the 2024 Noto Peninsula earthquake (Ministry of Land, Infrastructure, Transport and Tourism, Japan, 2014).To date, such tsunami scenario faults have been proposed to rupture individually, and triggering or simultaneous ruptures along with other adjacent scenario faults with considerably different fault geometries have not been specifically supposed.However, multiplex rupturing across a complex fault network or segmented multi-fault triggering has been reported for earthquakes in variable tectonic environments (e.g.Kuge et al., 1996;Satriano et al., 2012;Fan et al., 2017;Nissen et al., 2016;Hamling et al., 2017;Lay et al., 2018;Yamashita, Yagi, & Okuwaki, 2022;Vasyura-Bathke et al., 2024).This can result in diverse rupture behaviors, including the rupture of faults with conjugate dip directions and apparent backward rupture propagation.
As shown in the Results section, the E2 and E3 occur almost simultaneously, but they are separated by the preceding earthquake swarm region around the hypocenter.This complicates the apparent rupture directions of the 2024 Noto Peninsula earthquake.That is, the rupture first propagates west from the hypocenter, and then, it apparently goes both west and east of the hypocenter, resulting in the separated E2 and E3 with a change of the dip direction across the hypocenter.The regions of the preceding earthquake swarm and the 2023 M 6.5 source area around the 2024 hypocenter did not completely act as barriers to stop the 2024 ruptures.This may add critical inputs when considering earthquake hazards in a context of short-term assessment after a series of notable events (e.g.geodetically detected aseismic deformation, intense earthquake swarms, and a major earthquake).The 2024 Noto Peninsula earthquake suggests the triggering of segmented fault ruptures with opposite dip directions.It is therefore important to comprehensively consider interaction with preceding seismic or geodetic activities when assessing future earthquake hazards.

Conclusions
We found a series of discrete rupture episodes on differently oriented fault segments for the 2024 Noto Peninsula earthquake.The 10-s-long quiet initial rupture domain overlaps the preceding earthquake swarm region, which hosted the irregular nucleation process that leads to the dynamic instabilities of the major ruptures that followed.The multi-scale rupture growth was controlled by the complex fault network and heterogeneous source environment that is highlighted by swarm activity and associated aseismic deformation.A possibility of co-ruptures of the segmented faults or triggering of isolated fault ruptures with different fault orientations can add critical inputs into better assessments of future earthquake hazards.

Figure 1 .
Figure 1.Summary of the study region of the 2024 Noto Peninsula earthquake.(a) The star denotes the epicenter.The beachball shows the estimated total moment tensor obtained by the integration of all the potency-rate density tensors.The black circles show the 24-hour aftershocks (Japan Meteorological Agency, 2024).The cross markers show the past major events.The shaded areas show the regions of the 2020-2023 earthquake swarms (pink) and the 1-week aftershocks of the 2023 M 6.5 earthquake (gray) based on the Gaussian kernel-density estimates (> 5% of the normalized density).The red lines are the active faults (AIST, 2012).The topography is from SRT-MGL1 tiles (NASA JPL, 2013).(b) The regional map of the study region.The rectangle shows the area of Fig. 1a.(c) The station distribution (triangles) used for the analysis.The star denotes the epicenter.The dashed circles show the epicentral distances.

Figure 2 .
Figure 2. (a) The potency density tensor distribution, calculated by the integration of the potency-rate density tensors with respect to time.The star shows the hypocenter, which gives the initial rupture point.The rectangle outlines the regions of model faults.The thick black lines represent the top of the model faults.The black contour outlines potency every 1.1 m.(b) Potencyrate density evolution projected along the line of each model strike.The star represents the initial rupture point.The dashed rectangles highlight the notable rupture domains that are displayed in Figs. 3 and 4. The lines of V P and V S show the P and S wave speeds from the near-source structural model (TableS1).-12-

Figure 3 .
Figure 3. Selected snapshots of the strike orientations, extracted from the corresponding potency-rate density tensors.The bar represents the strike orientation, which is one of the two possible nodal planes of the best-fitting double couple solution that minimizes the inner product of fault-normal vectors of the candidate plane and the corresponding model plane.The bar length is scaled with the potency rate.The full snapshots of the potency-rate density tensor distribution can be seen in Movie S1.

Figure 4 .
Figure 4. (a) The evolution of strike orientations for the twin and single fault models, focusing on the E1-E3 rupture episodes.(b) The summary of potency-rate density distributions and the corresponding strike orientations.The colored contours show the areas of > 50% of the maximum potency rate from 15 to 36 s in the western segment and from 26 to 33 s in the eastern segment.The bar represents the strike orientation for the maximum potency rate in the corresponding time window.The inset shows the moment rate function.

Figure S1 .
Figure S1.Model geometries for the (a) twin and (b) single fault models.The blue and pink dots are the locations of source elements.The star shows the initial rupture point.The black circles show the 24-hour aftershocks (Japan Meteorological Agency, 2024).(c) The 24-hour aftershocks (Japan Meteorological Agency, 2024) colored by depth.

Figure S3 .
Figure S3.Comparison of potency density tensor distributions for the (a) twin and (b) single fault models.The figure style is the same as Fig. 2. The black contour outlines potency every 1.1 m for both the models.

Figure S4 .
Figure S4.The upper panels show the selected snapshots of the potency-rate density tensors and the corresponding strike orientation distributions.The bar represents the strike orientation, which is one of the two possible nodal planes of the best-fitting double couple solution that minimizes the inner product of fault-normal vectors of the candidate plane and the corresponding model plane.The lower panel shows the summary of the selected snapshots of the potency-rate density distribution and the strike orientations.The figure style and legends are the same as Fig. 4. The full snapshots can be seen in Movie S2.

Figure S5 .
Figure S5.Comparison of solutions built upon the (a) twin and (b) single model faults.The black contours highlight the loci of > 50% of maximum potency rate.The figure style is the same as Fig. 2.

Figure S6 .
Figure S6.Comparison of spatiotemporal distributions of potency rate density and the corresponding strike orientation from the twin (blue) and single (pink) fault models.The figure style and legends are the same as Figs. 4 and S4.

Figure S7 .
Figure S7.Waveform fits of the twin fault model.The black and red traces are the observed and synthetic waveforms.The station code and channel, the maximum amplitude of observed waveform, the station azimuth (φ), and the epicentral distance (∆) are shown on the left of each panel.The bottom map is an azimuthal equidistant projection of the station distribution (triangle).The star shows the epicenter.The dashed lines are the epicentral distances at 30 • and 90 • .

Figure S8 .
Figure S8.Waveform fits of the single fault model.The black and red traces are the observed and synthetic waveforms.The station code and channel, the maximum amplitude of observed waveform, the station azimuth (φ), and the epicentral distance (∆) are shown on the left of each panel.The bottom map is an azimuthal equidistant projection of the station distribution (triangle).The star shows the epicenter.The dashed lines are the epicentral distances at 30 • and 90 • .

Figure S10 .
Figure S10.Spatiotemporal distribution of the strike orientations.(a) The abscissa is a distance from the initial rupture point along each model strike, and the ordinate is a strike orientation (a remainder of dividing the strike angle by 180 • ).The strike orientation is extracted from the resultant potency density tensor, selected from the one of the two possible nodal planes of the best-fitting double couple solution that minimizes the inner product of fault-normal vectors of the candidate plane and the corresponding model plane.The star represents the initial rupture point.(b) The abscissa is a time from the origin time, and the ordinate is a strike orientation.The dashed rectangles highlight the notable rupture domains that are displayed in Figs. 3 and 4. TableS1).