Nearby fault interaction is a crucial yet often elusive process in seismic hazard analysis, influenced by regional tectonic stress and fault geometry. Understanding this interaction is particularly challenging in crustal thrust fault systems, where earthquakes on one fault can trigger events on other segments or sheets, sometimes even down to detachment zones. While interactions between major head-to-head, conjugate thrust faults are poorly understood, the 2022 Chihshang earthquake sequence in eastern Taiwan presents a unique opportunity for investigation. Eastern Taiwan experiences high tectonic shortening rates and frequent earthquakes, making it an ideal location to study this phenomenon. The Chihshang sequence included a Mw 6.5 foreshock followed by a Mw 7.0 mainshock, both occurring within the Longitudinal Valley, an active suture zone. This suture zone lies between the Luzon Arc and the accretionary wedge of the Taiwan Orogeny. The east-dipping, sinistral strike-slip Longitudinal Valley fault (LVF) is considered the main plate boundary fault, with a history of significant earthquakes. However, the activity of the west-dipping Central Range fault (CRF), located beneath the western side of the valley, has been debated. While a large earthquake may have occurred on the CRF in 1908, it has been largely inactive since then. The 2022 Chihshang earthquakes, located west of the Longitudinal Valley, demonstrate that the CRF also accommodates present-day plate convergence, forming an active head-to-head, conjugate fault system with the LVF. This study leverages seismic, geodetic, and field geological data to investigate the interaction between these two faults.

Previous research has documented various types of fault interactions in different tectonic settings. Studies have shown coseismic fault models from dense GPS networks (Ohta et al., 2008), complex multifault ruptures (Hamling et al., 2017), and the triggering of seismic events on different fault segments within the same thrust sheet (Zhang et al., 2020). Other research has highlighted the triggering of shallow minor structures (Béon et al., 2017; Yang et al., 2021). However, well-documented examples of interaction between two major head-to-head, conjugate thrust faults are rare. The San Fernando and Oak Ridge faults in Los Angeles offer one of the few such examples (Stein et al., 1994; Yeats & Huftile, 1995; Dolan et al., 1995). These faults produced the 1971 San Fernando and 1994 Northridge earthquakes, respectively, but their interaction before these events was poorly understood due to the low slip rates and long recurrence intervals of large ruptures (Dolan et al., 1995). The study of interseismic crustal deformation in the Taiwan plate boundary zone using GPS observations has also been conducted (Hsu et al., 2009), providing crucial background information for the current study. Previous work on the Longitudinal Valley fault has included studies on active fault creep variations (Lee et al., 2003), geomorphic analysis of the Central Range fault (Shyu et al., 2006), and interseismic coupling models (Thomas et al., 2014). These studies provided the groundwork for understanding the complex tectonic setting in which the 2022 Chihshang earthquakes occurred.

This study combined multiple datasets and analyses to characterize the fault interaction: 1. **Field Survey:** A post-mainshock field survey mapped surface ruptures along both the CRF and LVF, documenting offsets and fracture patterns. 2. **Optical Image Correlation:** Sentinel-2 satellite imagery underwent optical image correlation using the COSI-Corr method to estimate horizontal ground displacements associated with the foreshock and mainshock. The N-S displacement fields were processed with non-local mean filtering and resampled using a quad-tree scheme based on SNR and local variances. 3. **InSAR Analysis:** L-band ALOS-2 SAR data were processed to produce interferograms for both the foreshock and mainshock. The interferograms were unwrapped, converted to line-of-sight (LOS) displacements, and manually corrected for unwrapping errors. The results were validated against GNSS data and resampled based on coherence and local variances. Areas near fault traces were excluded due to low coherence. 4. **GNSS Measurements:** A dense GNSS array provided data to determine coseismic displacements for both the foreshock and mainshock. High-rate GNSS (HR-GNSS) data were processed with the GipsyX/RTGx software, incorporating atmospheric delay corrections and phase center modeling. Coseismic displacements were calculated by comparing average positions before and after the events. Postseismic displacements within 34 days after the mainshock were modeled using a specific function to capture the rapid postseismic deformation. 5. **Coseismic and Postseismic Slip Inversion:** A co-seismic slip inversion model was created using InSAR and GNSS displacements. The model included both CRF and LVF. The foreshock was modeled with CRF alone, while mainshock modeling included both faults. A non-negative least squares method penalized dextral slip and normal faulting, and the slip distribution was regularized with a second-order Laplacian operator to avoid overfitting. A similar approach was used for 34-day postseismic slip inversion with GNSS data. 6. **Kinematic Rupture Process Inversion:** HR-GNSS data were used in a joint source inversion with GNSS coseismic static displacements to analyze the kinematic rupture process. A spectral-element method was used to calculate Green's functions, and non-negative least-squares were applied to the joint inversion. Various constraints including maximum rupture speed, damping, and smoothing were employed. 7. **Coulomb Stress Estimation:** Coulomb stress changes after the 1951 and 2022 earthquake sequences were calculated using slip models. Stress changes were calculated with kernels that relate slip to stress changes, considering rake angles and effective friction coefficients. The 1951 earthquake slip model, and the 2022 foreshock and mainshock slip models from this study, were used as stress change sources. Conservative values were used in the moment budget calculations. 8. **Moment Budget Estimation:** A first-order moment budget for the CRF and LVF was estimated over the past 120 years. Moment accumulation was calculated using long-term slip rates, shear modulus, fault area, and time. Moment release was calculated from major earthquake seismic moments and aseismic moment release from the shallow LVF creeping zone. The moment budget was set to zero at previous large earthquakes (1908 for CRF and 1951 for LVF). The impact of the stress shadow effect from both the 1951 and 2022 earthquakes were considered in the budget calculation.

The 2022 Chihshang earthquake sequence revealed significant insights into the behavior of the CRF and LVF: 1. **Surface Rupture and Displacement:** Surface ruptures were observed along both faults. The CRF rupture was predominantly to the north, while the LVF showed minor fractures with left-lateral motions and shortening. Optical image correlation confirmed the north-south offset diminishing southwards and eastward shift of the offset discontinuity toward the LVF. 2. **Coseismic Slip:** The foreshock slip mainly occurred on the CRF, south of the hypocenter. The mainshock involved significant slip on the CRF (up to ~2.4 m) with a bilateral rupture process, accompanied by shallow slip (<5 km) on the LVF near Chihshang. HR-GNSS data revealed southward propagation during the foreshock and northward propagation during the mainshock, along with shallow LVF rupture between Chihshang and Yuli. 3. **Slip Behaviors:** The CRF showed subtle or no afterslip at the up-dip side, contrasting with the significant afterslip on the LVF near Chihshang. This difference may be related to lithological variations. Deep postseismic slip on the CRF may represent viscoelastic flow, while shallow LVF slip suggests a frictionally transitional zone capable of slip during various earthquake cycle stages. The shallow part of LVF showed intermediate interseismic coupling, possibly affected by CRF. 4. **Fault Interaction:** The close proximity of CRF and LVF led to intense interaction. Moment release on the LVF during the mainshock (9–15% of total) started slightly later than CRF, indicating dynamic triggering by seismic waves. The long duration of LVF slip suggests potential influence from the second asperity on the CRF. Historical records suggest simultaneous rupture of both faults during major earthquakes (M>7), alternating as primary and secondary structures. 5. **Coulomb Stress Changes:** The foreshock increased Coulomb stress on the CRF, triggering the mainshock. The mainshock caused a significant stress drop on the LVF, suggesting potential inhibition of seismicity. The 1951 earthquake similarly induced stress drops on the CRF, resulting in a period of quiescence. 6. **Historical Earthquake Records:** Historical data (since 1900) show an out-of-phase pattern of moment release between CRF and LVF, with quiescence periods following major events on the opposite fault. This out-of-phase pattern is the first time reported in such a detailed way and time-scale. 7. **Moment Budgets:** Analysis of moment accumulation (using long-term slip rates) and release suggests that the 2022 earthquake released most of the accumulated moment on the CRF since 1908. However, a significant moment remains on the LVF since 1951, indicating a potential for another major earthquake in the future. The stress shadow effect on the moment budgets was also considered.

The findings of this study significantly advance our understanding of fault interactions in complex tectonic settings. The 2022 Chihshang earthquake sequence highlights the active role of both the CRF and LVF in accommodating plate convergence, contrary to previous assumptions prioritizing the LVF. The observed profound fault interaction, characterized by periods of quiescence and out-of-phase moment release, underscores the necessity for reassessing seismic hazards in the region. This out-of-phase behavior, previously reported over larger scales in other regions, is here documented in detail for the first time in a centennial time scale along two major faults. The observed stress shadow effect requires consideration in seismic hazard assessment, as it may prolong quiescence periods and influence the timing of future large earthquakes. The study's moment budget analysis provides a valuable first-order estimation of potential future seismic activity, highlighting the significant remaining moment on the LVF and the potential for another major earthquake in the coming decades. This research should encourage further detailed studies incorporating sophisticated models accounting for fault geometry, slip characteristics, and stress interactions for more accurate seismic hazard estimations.

The 2022 Chihshang earthquake sequence dramatically demonstrated the presence and significance of the head-to-head, conjugate fault system formed by the CRF and LVF. The profound fault interaction observed highlights the need for a revised understanding of plate convergence in eastern Taiwan and necessitates integrating nearby fault interactions into seismic hazard assessments. Future studies should focus on more detailed modeling to improve estimates of seismic risks by incorporating geometric complexity, slip characteristics, and stress interactions.

The study's moment budget estimations are based on a simplified model using first-order approximations and long-term average slip rates. The uncertainties in these values affect the precision of the budget calculations. Furthermore, the sparse GNSS stations in the Central Range and the limited resolution of geodetic data pose challenges in determining the detailed geometry of the CRF. While the study incorporated various datasets and methods, the complexity of fault interactions might not be entirely captured by the current models. Further research is needed to refine our understanding of the interaction dynamics.