Nearby Fault Interaction within the Double-Vergence Suture in Eastern Taiwan during the 2022 Chihshang Earthquake Sequence

11 Faults interact through stress changes induced by earthquakes and viscoelastic flow acting in the 12 lithosphere, but the process can be geometry-dependent and temporally-variant. The Chihshang


Introduction
Nearby fault interaction, a process critical for seismic risk analysis, depends on the regional tectonic stress fields and the geometric setting of the faults.In crustal thrust fault systems, an earthquake on a fault may induce seismic events on a different fault segment within the same thrust sheet [1][2][3] , stepping over to another thrust sheet [3][4][5] , and/or down to the low-angle detachment 3 .
Ruptures on shallow minor structures may also be triggered at the same time 5,6 .Of all geometric settings, the interaction between two major head-to-head conjugate thrust faults has the least known cases-by far the only well-documented example is probably the south-vergent San Fernando fault and the north-vergent Oak Ridge fault in the Los Angeles region, which produced the 1971 M7.1 San Fernando earthquake and the 1994 M6. 8 Northridge earthquake, respectively [7][8][9][10][11][12] .The total length of the double-vergence fault system may extend for more than 90 km by counting in two more south-vergent faults further to the west (the Santa Susana fault and the San Cayetano fault), and the closest distance between faults of opposite vergence is less than 3 km.Before the 1971 earthquake, there was no record of large earthquakes along this fault system except for a probable one in December 1812 12,13 .Given the low slip rates of faults and the long recurrence interval for large rupture events [12][13][14] , it is not certain when the next large earthquake may happen and what the rupture scenario may be.However, a recent earthquake sequence in eastern Taiwan, along with high tectonic shortening rates and frequent earthquakes, may offer insights regarding the nearby fault interaction in a head-to-head conjugate fault system.
The September 2022 Chihshang earthquake sequence includes two major shallow earthquakes in eastern Taiwan.The first one is a moment magnitude (Mw) 6.5 foreshock on September 17, followed by an Mw 7.0 mainshock 7 km to the north and 17 hours later on September 18 (moment magnitudes from autoBATS CMT at https://bats.earth.sinica.edu.tw/).These two earthquakes occurred within the narrow and elongated Longitudinal Valley, an active suture zone between the Luzon Arc on the Philippine Sea Plate to the east and the accretionary wedge of the Taiwan Orogeny along the Eurasian continental margin to the west (Fig. 1a).Global Navigation Satellite System (GNSS) observations show that around 30 out of the total 90 mm/yr northwestward plate convergence is absorbed across the valley 15 .Within the suture, the east-dipping, sinistral strike-slip Longitudinal Valley fault (LVF) is considered the main plate boundary structure located along the eastern edge of the valley 16 (Fig. 1b).Several Mw 6.7-7.3 earthquakes have been documented along this fault during the past one hundred years 17 .On the other side of the valley, a west-dipping Central Range fault (CRF) had been proposed but faced many doubts and skepticism about its activity due to its close distance to the LVF [18][19][20][21] (Fig. 1b).Although a magnitude ~7 earthquake may have possibly occurred on this fault in 1908 17,22,23 , the south-central segment of the CRF (south of 23.5°N) shows no sign of creeping [24][25][26] and has remained seismically inactive, compared with the LVF, for nearly a century.The quiescence is eventually ended by the 2022 earthquakes.
The 2022 earthquakes are located to the west of the Longitudinal Valley, with the mainshock close to the town of Chihshang (Fig. 1d).This event demonstrates that the CRF accommodates the present-day plate boundary shortening, and that the LVF and CRF together form an active head-to-head conjugate thrust system (Fig. 1b), similar to the San Fernando fault and the Oak Ridge fault in southern California.With this unique opportunity, we look into the interaction between the two faults during and after the earthquakes.We further assess the first-order moment budgets on these two faults by using seismic, geodetic, and geologic records.

Surface rupture and displacement field
Field survey immediately after the mainshock shows that most of the surface ruptures occurred 15 km north of the epicenter, close to the town of Yuli (Fig. 1d, 2).The rupture traces are 1 to 2 km east of the previously mapped CRF along the Yuli fault, which was first identified by minor surface ruptures during the 1951 Mw 7 earthquake 17,29,30 (Fig. 1c).The Yuli fault used to be considered as an en echelon stepover of the LVF 29 , but now we reattribute it to the larger CRF system (with "system" omitted hereafter for simplicity).Significant left-lateral and vertical offsets (west side up) are observed in this section (Fig. 2b, c), whereas to the south near the epicenter, field survey shows only diffuse deformation without noticeable offsets.On the LVF, minor fractures are distributed all the way from north to south showing transpressional left-lateral motions (Fig. 2a, d,    e).These field ruptures are comparable to the results from optical image pixel offset tracking using Sentinel-2 images: the sharp 2-m N-S offset in the north diminishes southwards, and the offset discontinuity moves eastwards towards the surface trace of the LVF (Fig. 1e).The topographic relief in the mountain ranges on both sides of the valley limits the resolvable pixel offsets away from the fault, so interferograms from synthetic aperture radar images acquired by the L-band Advanced Land Observing Satellite-2 (ALOS-2) provide additional constraints for the coseismic deformation (Fig. 1c, d).The ground displacement during the foreshock displays an anti-symmetric pattern across the fault.During the mainshock, the ground west of the CRF has moved ~1 m towards the satellite, while the ground to the east has moved ~0.1 m away from the satellite.This asymmetric pattern indicates significant reverse motion along a west-dipping fault plane, compatible with the field observations and focal mechanisms (Fig. 1d; solutions from the Real-Time Moment Tensor Monitoring System 31 at https://rmt.earth.sinica.edu.tw/).The foreshock and the mainshock are each located at the northern and southern edge of the corresponding deformation patches, suggesting bilateral rupturing processes that occur sequentially with a 17-hr time delay.
In the meanwhile, more than 30 GNSS stations operate continuously along the Longitudinal Valley.During the foreshock, a dense GNSS array near Chihshang displays a sharp change of N-S displacements from the hanging wall to the footwall of the CRF (Fig. 1c).During the mainshock, the coseismic horizontal and vertical displacements west of the CRF reach 0.65 m and 0.95 m, respectively (Fig. 1d, e).East of the CRF, the maximum horizontal displacement also reaches 0.72 m, with the maximum ground subsidence of 0.37 m.Along the dense GNSS array near Chihshang, station TAPE and TAPO, which are ~700 m apart across the LVF, show opposite horizontal motions during the mainshock (Fig. 2f, g).Such a short-wavelength deformation, together with surface ruptures observed along the LVF (Fig. 2a, d, e), suggests that coseismic slip occurred in the shallow part of the LVF during the mainshock.mainshock.The northward motion at TAPO suggests that the station is located at the footwall of the CRF.However, the vertical motion at TAPO is up, different from the subsidence pattern observed in all other stations east of the LVF (same direction as CE2A).This anomalous motion suggests coseismic shallow slip on the LVF.

Coseismic slip and rupture process
Given that surface ruptures are observed simultaneously along the CRF and LVF faults, the fault model is designed with a west-dipping plane for the CRF, with a uniform dip angle of 60° according to the focal mechanism, and an east-dipping fault plane, with its listric geometry and surface trace following the parameters for the LVF in the Taiwan Earthquake Model 28 .In the northern section near Yuli, the surface trace of the CRF follows the observed ruptures and the discontinuity in pixel offset tracking results (Fig. 1e).In the southern section, to account for the lack of surface ruptures on the CRF and the eastward shift of the discontinuity in the N-S displacement field, the up-dip end of the CRF is shifted eastwards and submerges under the LVF.
Our model result shows that the CRF alone is sufficient to explain the coseismic displacements during the foreshock.In contrast, both the CRF and LVF are required to model the coseismic displacements during the mainshock, especially the near-field deformation across the LVF observed by HR-GNSS.( We use satellite data and GNSS displacements to constrain the coseismic slip distribution on the two faults (Fig. 3a).During the foreshock, most of the slip occurs to the south of the foreshock hypocenter.During the mainshock, the rupture forms two major asperities: the first one is right to the north of the epicenter at a depth of 5-15 km; the second asperity starts further north from a greater depth, and extends to the surface north of Yuli.The slip peaks at ~2.4 m, with an average rake of 53°.The total slip from the foreshock to the mainshock forms a contiguous patch over 60 km long on the CRF.On the LVF, minor and predominantly left-lateral slip takes place mostly at the shallow depth (<5 km) near Chihshang during the mainshock.Our model also indicates another shallow slip patch on the LVF near Yuli.Although there is no direct evidence from nearby GNSS stations like those near Chihshang, two bridges near Yuli were damaged due to substantial surface offset along the LVF (Fig. 2a).
To better understand the rupture process, 1-Hz sampling high-rate GNSS (HR-GNSS) position time series are directly used to invert for the kinematic rupture with the constraint of final displacements equal to the GNSS coseismic static displacements.The result shows that the rupture propagates on the CRF predominantly southwards during the foreshock, and it terminates at ~5 km depth (Fig. 3c, d).The lack of shallow slip is slightly different from the model constrained with satellite data (Fig. 3a), suggesting possible shallow afterslip on the CRF after the foreshock.
Actually, according to the GNSS position time series, creeping continued for some time after the foreshock near Chihshang (fig.S6).A few hours later, the mainshock rupture initiated close to where the foreshock rupture terminated and propagated unilaterally northwards, accompanied by the shallow ruptures on the LVF between Chihshang and Yuli (Fig. 3e, g).The seismic moment release estimated from the two aforementioned modeling strategies is 5.7-9.7×10 18Nꞏm (equivalent to Mw 6.4-6.6) for the foreshock and 3.3-3.6×10 19Nꞏm (equivalent to Mw 7.0) for the mainshock.

Slip behaviors on the CRF and LVF
This earthquake gives us a glimpse into the slip behaviors on the long-debated CRF.Based on geodetic records, the shallow part of the CRF displays no clear sign of interseismic creeping from Yuli to Chihshang 25 and was ruptured during the earthquake.To study the variation of slip behaviors along the down-dip direction, we model the GNSS postseismic displacements within the first 34 days following the earthquake as afterslip to obtain the approximate spatial distribution (see fig.S16 and table S3 for scenario tests between single and dual faults).The result shows subtle to none afterslip at the up-dip side of the CRF rupture zone despite the apparent coseismic shallow slip deficit (Fig. 4).However, we caution that the lack of significant shallow afterslip may as well result from sparse GNSS distribution close to the CRF north of Chihshang.On the other hand, notable shallow afterslip occurs on the LVF near Chihshang.Such a contrast is likely attributable to the lithological difference between the younger, non-metamorphic mélange east of the LVF and the older, metamorphic mélange west of the CRF 32 .The down-dip side of the CRF rupture zone is bounded by afterslip of ~0.8 m extending to a depth of 30 km (Fig. 4).The deep postseismic slip is mostly located in the aseismic zone beneath the southeastern Central Range, which is believed to be undergoing ductile deformation at low viscosity due to high heat flows and partial melting 33,34 .
Only a few aftershocks were recorded in this region within two months of the earthquake, suggesting that viscoelastic flow may play an essential role in postseismic deformation in the vicinity of the fault and from this depth downwards.Therefore, part of the inverted afterslip may be a proxy of the viscoelastic flow excited by the deviatoric stress changes after the earthquakes.
Note that the southern end of the afterslip patch seems to have a larger amount of slip, which is possibly a result of uneven station distribution.Incorporating satellite-based observations may help verify this anomaly.
On the LVF, the shallow part has been creeping constantly at high rates 19,35,36 and experienced coseismically-induced slip during the mainshock.Postseismic GNSS observations and modeling suggest that afterslip also occurs close to the surface near Chihshang (Fig. 4), despite a Coulomb stress drop of ~1 MPa after the earthquake (Fig. 5c).The fact that slip can occur on the shallow part of the LVF during different stages of the earthquake cycle is rather unique and perplexing.It implies that more strain is released during and after the earthquake, even though the strain is being released constantly in the interseismic period.This phenomenon suggests that the shallow part of the LVF is probably a frictionally transition zone, which is a conditionally stable region 37 that can also promote the propagation of seismic ruptures depending on the stress and physical conditions 38 .Studies have shown strong ground shaking during large events may potentially change the elastic and frictional properties in the fault zone 39 .Under what exact conditions the shallow part can undergo all-time slip requires further explorations from laboratory experiments and numerical simulations.
A previous geodetically-determined interseismic coupling (ISC) model for the LVF indicates intermediate coupling at the shallow part and low coupling at the deeper part near Yuli 36 (Fig. 4).The spatial extent and orientation (N20˚) of the shallow coupling patch agree well with those of the 2022 coseismic slip patch, suggesting that the ISC model has indeed captured the elastic strain associated with fault coupling, albeit the main source is the CRF.The actual coupling for the LVF, from shallow to deep, is therefore in question given the nearby influence from the CRF.In addition to geometric unknowns (single vs. dual faults as well as fault intersections), the short geodetic baselines in the Longitudinal Valley may limit the resolution at depths.An updated ISC model with both faults considered and offshore geodetic constraints will be needed to obtain realistic estimates of the coupling status on the faults.Gray circles are aftershocks between the foreshock and mainshock.(c) Coulomb stress changes due to the slip in the mainshock.Gray circles are the aftershocks after the mainshock.A rake of 45° is adopted for the CRF and LVF as a receiver fault.The friction is assumed to be 0.4 40 .Note the significant stress decrease on one fault due to the slip on the other.Earthquake catalog is from Taiwan GDMS.

Nearby fault interaction during earthquakes
The close distance between the two head-to-head conjugate thrust faults leads to intense and complicated interactions during seismic events.The moment release on the LVF is 3.2-4.9×10 18 in the mainshock, accounting for 9-15% of the total moment release.According to the source time functions for the mainshock (Fig. 3b), the moment release rate on the CRF rose quickly, peaked at 9 s, and decreased stepwise till 30 s.The moment release on the LVF started 2-3 s later, increased gradually all the way till 20 s, followed by a sharp drop at 20-23 s and another tiny peak at 25-26 s.
The slightly-delayed onset of moment release on the LVF suggests that the shallow slip was dynamically triggered by seismic waves, most likely by the strong S wave or Rayleigh wave 38 .The long duration for the induced slip on the LVF is particularly intriguing.The slip even picked up more momentum during the later phase of the mainshock between 15-20 s while the moment release on the CRF has already decreased from the peak.One possible explanation is a larger excitation when the S wave or the surface wave from the second asperity traveled northeastwards across the LVF.In view of dual ruptures, the same phenomenon has actually been recorded previously during the 1951 earthquake, whereas the LVF hosted the main slip (Fig. 1c).The records from both the 1951 and 2022 earthquakes suggest that the CRF and LVF are likely to rupture simultaneously during major earthquakes (Mw>7), with the faults alternating between the primary and the secondary structures.Such a tendency needs to be considered when evaluating the energy budget of seismic events.

Static stress changes in a closely-spaced fault system
As what would be expected from the perspective of typical Coulomb stress transfer, the ruptures on the CRF during the foreshock caused a Coulomb stress increase by 0.5 MPa to the northern edge of the slip patch, which then led to a sequence of aftershocks and eventually the mainshock (Fig. 5b).After the mainshock, a ring of increased Coulomb stress circled the coseismic slip patch on the CRF, with the stress increase up to 2 MPa at the up-dip and down-dip ends.On the LVF, however, a continuous patch with stress reduction extends through the entire fault section at 2-20 km depth, with the maximum decrease of 1 MPa near the hypocenter (Fig. 5c).This large stress drop suggests that the seismic potential on the LVF may be impeded, at least in principle, subsequent to the 2022 earthquakes.
Similarly, the 1951 coseismic ruptures 16 induced a Coulomb stress drop of 1-2 MPa on the CRF (Fig. 5a).Although these values may be overestimated given the uncorrected postseismic effects, the pattern and sign of stress changes can still be instructive.They correspondingly suggest that the seismicity on the CRF may have been suppressed following the 1951 earthquakes.
The temporary shutdown of seismicity caused by slip on the other fault in a head-to-head conjugate thrust fault system can also be found in historical earthquake records.When viewed from the entire length of 140 km (Fig. 6a), historical earthquakes since 1900 with magnitude larger than 6 and depth shallower than 30 km 17,23 show a period of quiescence following the occurrence of a major event on the other fault-a 37-year-long quiescence on the LVF after the 1908 Mw 7 earthquake on the CRF, and a 62-year-long quiescence on the CRF after the 1951 earthquake sequence on the LVF (Fig. 6b).A larger moment release on one fault likely results in a longer quiescence period on the other fault.This relationship agrees with the sizes of stress drop estimated after the 1951 and the 2022 earthquakes (Fig. 5).It also suggests that for the head-to-head conjugate thrust faults with lower slip rates, such as those in the Los Angeles region (3-7 mm/yr for the Cayetano fault 12,13 and 3-12 mm/yr for the Oak Ridge fault 10,12 as compared to 29-38 mm/yr for the LVF and 6-11 mm/yr for the CRF; see table S4), the stress drop imparted by an Mw>7 earthquake on the nearby fault may lead to an even longer quiescence.This effect may take part in the centurieslong seismic lull in the Los Angeles area, especially in between the short-range thrust faults 12 .
Arguably, as a result of strong stress shadow effect, the moment release forms a concaveup curve on the CRF and a concave-down curve on the LVF (Fig. 6c, d).The two curves together indicate that the moment releases on the two faults are out of phase.Such an out-of-phase pattern has also been reported over a larger spatiotemporal scale across different fault systems in Southern California 41 and in central Italy 42 .It is however the first time that the out-of-phase seismic bursts are reported along two major faults in a hundred-year time scale.In the long run, viscous flows in the lower crust may weigh in to add more complexity to the interaction between the faults and with other structural systems in the far range 42 .Further physics-based modeling will be needed to examine the stochastics and stationarity of the out-of-phase pattern.

Moment budgets and seismic risks
A large earthquake like the 2022 mainshock raises the question of how many more large events to be expected along the CRF and LVF.While a detailed interseismic coupling model may be preferred to answer such a question, such modeling has a large degree of freedom in this area: sparse GNSS stations in the Central Range, poorly-constrained geometry of the CRF, and the limited spatial coverage of geodetic data around the Longitudinal Valley.Alternatively, we invoke centennial earthquake records and millennial geologic rates to estimate a first-order moment budget for the two faults over the past 120 years.The moment budget consists of moment accumulation and moment release.To estimate the moment accumulation, a 140-km × 30-km plane is assigned to represent the full length of each fault.The long-term slip rate  of 29.0-38.0mm/yr and 6.7-11.0mm/yr is assigned to the LVF 43,44 and the CRF 20,28 , respectively, based on the long-term uplift of fluvial terraces on each side (table S4).The accumulated moment can be estimated by using  ∆, where  is shear modulus (30 GPa),  is the fault area, and ∆ is the year since 1900.The moment release can be estimated by summing up the seismic moments from major earthquakes (Fig. 6a) and the aseismic moment release from the shallow creeping zone on the LVF (table S4).For simplicity, the moment budgets are set to zero at the previous large earthquakes, i.e., the 1951 earthquake sequence for the LVF and the 1908 earthquake sequence for the CRF (Fig. 7b).The result shows that the 2022 earthquake sequence seems to release most of the moments accumulated on the CRF since 1908.The estimated current budget is around −0.7~5.9×10 19Nꞏm given the uncertainties in the long-term slip rate.The surplus can be accommodated by small historical events (M<6) that are not accounted in the estimate, the postseismic transients, and/or the future earthquakes that fill up spatial gaps (Fig. 7a).Note that the spatial gap north of 23.75°N is surrounded by small rupture patches associated with Mw 6.1-6.4 earthquakes, and hence the gap may serve as a possible candidate to host future large earthquakes.
On the other hand, a considerable amount of moment has continued to accumulate on the LVF since 1951.The current budget, with the moment release from shallow interseismic creeping subtracted, is already 2-3 times the moment released by 2022 earthquakes (Fig. 7b).The current surplus reaches 1.7-3.1×10 20Nꞏm, equivalent to an earthquake of Mw 7.4-7.6.On the map view, Mw 6-6.8 earthquakes have occurred sporadically along the deeper part of the LVF, while the shallow part (<10-15 km) has not experienced major seismic events since 1951 (Fig. 7a).Although smaller earthquakes and interseismic creeping may consume part of the accumulated moment, another series of major earthquakes with the same amount of moment release as the 1951 earthquake sequence is likely to take place within 45-84 years given the current budget accumulation trend.
The influence of stress shadow on the CRF and LVF imposed by the 1951 and 2022 earthquakes, respectively, is also included in the moment budget (see "Methods").The stress shadow effect on the LVF's budget curve is difficult to visualize due to another Mw 6.8 earthquake in 2022, while the effect is visible on the CRF's budget curve as a small offset (Fig. 7b).This offset indicates a drop of stress level on the CRF and seismicity shutdown afterwards, as revealed in the post-1951 quiescence.Although this is a rough estimate, we show that nearby fault interaction may modulate the moment budget in a head-to-head conjugate thrust fault system such as the CRF and LVF.Fig. 7. Rupture patches and moment budget for the CRF and LVF.(a) Map view of the rupture patches for Mw>6 earthquakes [45][46][47][48][49][50] .The northern extension of the CRF is based on the distribution of background seismicity 27 and rupture patches 50 .The up-dip end to the north and to the south is likely blind, overridden by the LVF.(b) Moment budget of the CRF and LVF.The uncertainty estimation comes from the uncertainties in long-term slip rates.

Summary
The presence of a head-to-head conjugate thrust fault system along the Longitudinal Valley suture was not completely obscure given the historical earthquakes and the geomorphic evidence of the CRF.However, the system did not draw much attention in the past, partly due to the high activity exhibited by the LVF.The 2022 Chihshang earthquake sequence demonstrates for the first time that either fault in the system is capable of generating earthquakes of Mw 7 or greater.Profound fault interaction during and after the earthquake has caused periods of quiescence and an out-ofphase temporal pattern of moment release between the two faults.This result not only urges us to revise the conventional thinking of plate convergence across eastern Taiwan, but also signals the importance of incorporating nearby fault interaction when assessing seismic hazards in similar tectonic settings.Future modeling needs to account for the geometric complexity, slip characteristics, and stress interactions between the faults in order to obtain realistic estimations.

Optical image correlation of Sentinel-2 image
Sentinel-2 level 2A (L2A) data (https://doi.org/10.5270/S2_-znk9xsj)acquired between 2022/08/23 and 202209/22 are used to estimate the horizontal ground displacements associated with both the foreshock and the mainshock.The L2A product includes orthorectified blue (B2), green (B3) and red (B4) bands at 10-m spatial resolution.To derive the ground displacement fields, the pre-and post-earthquake images were each converted to a gray-scale image using the linear combination of B2, B3 and B4 bands.The Co-registration of the Optically Sensed Images and Correlation (COSI-Corr) method 51 was used to retrieve the sub-pixel displacements between the two images.We utilized two window sizes, 64-km × 64-km pixels and 32-km × 32-km pixels, for the correlation, and set the sliding step to 1 pixel.The corresponding resolution for the displacement fields is 10 m.We set the robustness iteration to 4 and the frequency mask threshold to 0.98.The output of COSI-Corr includes three map layers: E-W, N-S displacement field and the signal-tonoise ratio (SNR) layer (fig.S2).To further decrease the noise in E-W and N-S displacement fields, we apply the non-local mean filter with the noise parameter H set to 1.6.The final results are validated with the surface offsets collected from the field survey.Because of the high noise level and the small displacement values, we discard the E-W displacement field.The N-S offsets are then resampled using a quad-tree scheme based on the data SNR values and the local variances.A subset of ~2250 points are selected for the kinematic modeling.

ALOS-2 L-Band SAR interferometric data
L-band ALOS-2 PALSAR-2 SAR data in ScanSAR mode acquired on 2022/04/03 and 2022/09/18 from descending path 27 are used to produce the interferogram for the foreshock.The post-event image was acquired 6 hours after the foreshock, so some post-seismic displacements are included in the image (fig.S3).ALOS-2 SAR data in Stripmap Fine mode acquired on 2022/05/19 and 2022/09/22 from ascending path 137 are used to produce the interferogram for the mainshock.
The post-event image was acquired 82 hours after the mainshock.A different track, ascending path 139 was available ~10 hours after the mainshock, but the images suffer strong ionospheric disturbance and were not able to provide as clear fringe patterns as those in path 137.The interferograms are produced using the standard Interferometric synthetic aperture radar Scientific Computing Environment package (ISCE2, https://github.com/isce-framework/isce2), w i t h t h e embedded ionospheric correction module for L-band wide-swath SAR data 52 .The unwrapped interferograms are converted to line-of-sight (LOS) displacements.Because of low coherence in the mountain areas, multiple unwrapping errors are observed in the interferogram that covers the mainshock.We manually correct for the unwrapping errors and validate the results with the coseismic displacements from GNSS observations projected to the LOS directions (fig.S4).The result shows a better agreement between the InSAR and the GNSS LOS displacements.The two interferograms are then resampled in a non-uniform scheme based on the coherence value and the local variances.For the foreshock interferogram, we crop out the region contaminated by ionospheric noises and resample the data to 1028 points following a quad-tree scheme based on the data coherence and local variances.The mainshock interferogram is resampled to 5797 points using the same method.

Processing of high-rate GNSS data
Given that the two earthquakes were less than one day apart, in order to precisely determine the coseismic displacements, high-rate GNSS observations from multiple hosting organizations were collected to estimate the instantaneous positioning, including the Institute of Earth Sciences, Academia Sinica (IES), Central Weather Bureau (CWB) and National Land Surveying and Mapping Center from Ministry of Interior (MOI).The GipsyX/RTGx (Real Time GIPSY) software package 53 was used to estimate the coordinates in precise point positioning (PPP) method with the precise final, non-fiducial daily products of orbit positions and clock offsets from the Jet Propulsion Laboratory archives.The Vienna Mapping Function 54 and antenna calibration provided by the National Geodetic Survey in National Oceanic and Atmospheric Administration were used to reduce atmosphere delays and to carry out phase center modeling.
The coseismic displacement field for the foreshock was estimated from the position difference between the averages of 1-hz solutions 30 seconds prior to and 30-60 seconds posterior to the earthquake.The maximum horizontal displacement is measured in the station ERPN, with a southwestward horizontal displacement of ~0.19 m.The coseismic displacement field of the mainshock was estimated from the position difference between the averages of 1-hz solutions one minute prior to and 1-2 minutes posterior to the mainshock.The maximum horizontal displacement is measured in the station Yul1, with a SSW horizontal displacement of ~0.75 m and a vertical displacement of almost one meter.The locations and time series of all the high-rate GNSS sites used in the modeling are shown in fig.S5 and S6.
For the postseismic displacements, we processed the position GNSS data within 34 days after the mainshock.The postseismic position time series   are fitted by the following function   34,55 : where A is the amplitude of postseismic transient,  is the time since mainshock,  is a dimensionless parameter related to the velocity-strengthening behavior of faults in response to stress changes, and  is a characteristic time scale of transient deformation.Here we use  = 5.5 and  = 0.1 years, determined by a grid search, to capture the rapid postseismic deformation.

Co-seismic and postseismic slip inversion
We include both the CRF and LVF in our model to image the coseismic slip (fig.S7).In the inversion, we invert the foreshock coseismic displacements with the CRF solely, while the LVF is included to jointly explain the mainshock observations.For the CRF, the fault trace near Yuli is illustrated by the coseismic surface rupture along the Yuli fault (a branch of the CRF).The surface rupture of the Yuli Fault vanishes southwards due to the overriding of the LVF.Therefore, we take advantage of the planar fault searched by the coseismic displacements of the foreshock, including GNSS and InSAR, to extend the CRF southwards to 22.9°N.Further south, we adopt the CRF previously proposed in the Taiwan Earthquake Model (TEM) 28 .Given the constraints from regional CMT solutions of the Chihshang sequence 31 and the aftershocks of the 2006 Mw 6.1 earthquake 45 , we set the west-dipping angle of 60° for the CRF from 0-20 km depth in the model.On the other hand, the listric geometry of the LVF is well documented in the TEM 28  To accommodate geodetic data with various temporal spans, we design an inversion scheme to simultaneously invert the coseismic slip during the foreshock and mainshock.The matrix form with linear observation equations for multiple geodetic data is formulated as where each  is a matrix and each  and  are data and model vectors, respectively.The subscripts of  and  denote the data type (D27: ALOS-2 descending track 27; A137: ALOS-2 ascending track 137; S2: Sentinel-2; G1: GNSS of the foreshock; G2: GNSS of the mainshock). and  with superscripts "CRF" or "LVF" represent the slip on the corresponding fault and the corresponding green's function with a Poisson ratio of 0.25 56 , respectively. with subscripts "fore" and "main" describes the fault slip during the foreshock and mainshock.Note that we limit the slip to south of Chihshang during the foreshock and north of 23°N during the mainshock to confine the slip distribution.In addition, we simultaneously fit a bilinear spatial ramp in the InSAR data to assimilate long-wavelength orbital errors that cannot be fully eliminated from data preprocessing.
The superscript "R" on  and  stands for a ramp and the subscripts "D27" or "A137" refer to the corresponding InSAR image.The relative weight between the InSAR and GNSS data is scaled to fit each dataset without a significant decrease of root-mean-square error (RMSE) to a specific data (fig.S8), while the Sentinel-2 data is weighted one-tenth of the InSAR data due to the relatively large perturbations in the data (approximately 10 times of the InSAR data).
Considering the sinistral oblique characteristics of the two faults, we perform non-negative least squares to penalize dextral-slip and normal faulting in the inverse problem.To avoid overfitting the data, we regularize the slip distribution on faults by where  is a discrete second-order Laplacian operator.The slip on the CRF and LVF share a unified hyperparameter  .The hyperparameter  is determined by the trade-off curve between misfit and model roughness (fig.S9).We also penalize the sum of the slip during the foreshock and mainshock individually by where  and  in the matrix are row vectors full of 1 and 0, respectively, and the hyperparameter  is the weight of the penalization.This is a weak penalization to suppress the slip on the subfaults with low resolution.The  and  adopted in our inverse model is 0.5 and 0.008, respectively.The RMSE of each dataset is documented in table S1.
The same approach is adopted to invert for the 34-day afterslip on the CRF and LVF with GNSS data (fig.S15).To accommodate the broad postseismic displacements recorded by GNSS in western Taiwan, we extend the CRF down-dip to a depth of 30 km.It is worth mentioning that we do not penalize the early afterslip on the coseismic patches in our model.The resulting model provides a RMSE of 9 mm for GNSS data (fig.S16, table S3).Table S1.RMSE (in cm) for coseismic displacement data fitting in different model scenarios Table S2.RMSE (in cm) for HR-GNSS data fitting in different model scenarios Table S3.RMSE (in cm) for postseismic data fitting in different model scenarios Table S4.Parameters used in moment budget estimation

Fig. 1 .
Fig. 1.Surface displacements of the September 2022 Chihshang earthquake sequence.(a) Tectonic setting of the Longitudinal Valley suture zone.CeR: Central Range; CoR: Coastal Range.(b) Cross section view showing the distribution of background seismicity between 1990 and 2020 (light gray dots) 27 .Relocated aftershocks prior to November 17 are in black dots (catalog from Taiwan Geophysical Database Management System, Taiwan GDMS, https://gdmsn.cwb.gov.tw/).See (d) for profile location.(c) Coseismic interferogram and GNSS horizontal displacements for the foreshock.Traces of the Central Range fault and the Longitudinal Valley fault (LVF) are based on the Taiwan Earthquake Model 28 .Traces of the 1951 surface ruptures on the Yuli fault (YLF) and the LVF are based on previous field surveys 29,30 .LOS: radar line-of-sight direction.(d) Coseismic interferogram and GNSS horizontal displacements for the mainshock.(e) Cumulative N-S offsets from the foreshock to the mainshock.Vectors represent the vertical displacement for the mainshock.

Fig. 2 .
Fig. 2. Field photos and geodetic observations of the ruptures along the CRF and LVF.(a) Map showing locations of the surface ruptures and field photos along (b-c) the CRF and (d-e) the LVF.(f) Coseismic GNSS horizontal displacements for the mainshock near Chihshang.(g) North and vertical component of the HR-GNSS position time-series at station TAPE (30 m west of the LVF), TAPO (540 m east of the LVF) and CE2A (2.3 km east of the LVF).FS: foreshock; MS:

Fig. 3 .
Fig. 3. Coseismic slip during the September 2022 Chihshang earthquake sequence.(a) Overview of the kinematic model constrained from satellite imagery and GNSS coseismic displacements.Brown contours and colored patches represent the slip during the foreshock and mainshock, respectively.Gray circles are aftershocks between the foreshock and mainshock.YL: Yuli.CS: Chihshang.(b) Source time functions obtained from the rupture process modeling constrained by HR-GNSS data and GNSS coseismic displacements.Gray and black plots represent the foreshock and the mainshock on the CRF, respectively, whereas orange curve shows contributions of the LVF during the mainshock.Records from two of the HR-GNSS sites (KUA2 and FUDN) are shown in (a).Black and magenta curves are observed and modeled records.(c-d) Rupture process for the foreshock.Hatched areas indicate the fault plane excluded from the inversion.(e-g) Rupture process for the mainshock.

Fig. 4 .
Fig. 4. Coseismic slip (during the foreshock and mainshock) and early afterslip (34 days after the mainshock) of the 2022 Chihshang earthquake sequence.Contour interval is 0.4 m on the CRF and 0.1 m on the LVF.The fault plane has been extended to 30 km at depth for the afterslip inversion.The interseismic coupling (ISC) model was constructed by using InSAR and GNSS observations with the geometry of the LVF only 36 .Gray dots are relocated background seismicity (M>1.5) between 10 and 30 km, 1990-2020 27 .Thick black dots are aftershocks (M>2) between September 18 and November 17 from Taiwan GDMS.YL: Yuli.CS: Chihshang.

Fig. 5 .
Fig. 5. Coulomb stress change after the 1951 and 2022 earthquake sequence.(a) Coulomb stress change on the CRF imparted by the slip on the LVF (slip model from 16 ) during the 1951 earthquake sequence.(b) Coulomb stress changes due to the slip in the September 2022 foreshock.

Fig. 6 .
Fig. 6.Historical earthquake records in the Longitudinal Valley since 1900.(a) Spatial distribution of historical earthquakes since 1900 with M>6 and depth <30 km 17,23 .Red and blue stars mark earthquakes attributed to the CRF and LVF, respectively.(b) Spatial and temporal distribution of historical earthquakes and the two quiescence periods.(c) Moment release by individual events on the CRF.(d) Moment release by individual events on the LVF.
. The east-dipping angle gradually decreases from 75° at the surface to 45° at a depth of 20 km in our model.The fault geometry is discretized into parallelogram subfaults approximately 2-km by 2-km in size, yielding 492 and 480 patches for the CRF and LVF, respectively.