Strong asymmetry in near-fault ground velocity during an oblique strike-slip earthquake revealed by waveform particle motions and dynamic rupture simulations

The 2022 Mw 7.0 Chihshang earthquake (Taiwan), captured by near-fault strong-motion seis-mometers, high-rate GPS, and satellite imagery, offers a rare opportunity to examine dynamic fault rupture in detail. Using dynamic rupture simulations, we investigate the particle motions recorded at near-fault strong-motion and 1 Hz GPS stations surrounding the main asperity. Some of these stations were as close as 250 m from the fault trace as determined by sub-pixel correlation of Sentinel-2 images. Our model reproduces the observed strong asymmetry in the ground motions on either side of the fault rupture, which results from along-dip spatial variability in rake angle on the steeply-dipping fault (70°) at shallow depth (2 km). Observed near-fault, pulse-like fault-parallel ground velocity larger than fault-normal velocity can be explained by a model with a sub-shear rupture speed, which may be due to shallow rupture propagation within low-velocity material and to free surface reflections. In addition, we estimate a slip-weakening distance D c of ~0.7-0.9 m from strong-motion seismogram recorded at Station F073, which is located ~250 m from the fault rupture, and the results of dynamic rupture modeling. The inferred D c is similar to other empirically derived estimates found for crustal earthquakes. These results have important implications for near-fault ground-motion hazard.


Introduction
Near-fault ground motion records of large magnitude earthquakes (Mw >7.0) typically contain strong faultparallel velocity pulses associated with permanent displacement of the ground surface.This feature is commonly referred to as 'fling step' in the earthquake engineering literature (Abrahamson, 2001;Kalkan and Kunnath, 2006).Observations adjacent to the surface fault rupture (<4 km) have been documented for large strikeslip earthquakes such as the 1999 Mw 7.3 Landers, 2002Mw. 7.9 Denali, 2016Mw 7.0 Kumamoto, and 2016 Mw 7.8 Kaikōura earthquakes (Hall et al., 1995;Kaneko and Goto, 2022;Kaneko et al., 2017;Ellsworth et al., 2004), and for thrust events such as 1999 Mw 7.6 Chi-Chi and 2008 Mw 7.9 Wenchuan earthquakes (Wen et al., 2010;Ji et al., 2003).Yet, for oblique-slip earthquakes that occur on steeply-dipping faults, near-fault observations are not as common.It is not clear how oblique slip may affect near-fault ground velocities, and to what degree across-strike symmetry (e.g. for strike-slip earth-quakes) or hangingwall effects (e.g. for thrust earthquakes) will dominate the near-fault ground motions in these settings.
The 2022 Mw 7.0 Chihshang earthquake struck the east coast of Taiwan on September 18 and was preceded 16 hours earlier by a Mw 6.6 foreshock (Lee et al., 2023;Yagi et al., 2023) (Fig. 1).Rapid moment tensor (Lee et al., 2013) and USGS finite fault solutions indicate both events likely occurred on north-east trending structures within the Longitudinal Valley, with focal mechanisms showing oblique-left-lateral slip consistent with the local seismotectonic setting (Fig. 1).The earthquake is consistent with the kinematics of oblique collision between the Luzon Volcanic Arc and Central Range in Taiwan, which gives rise to both thrust and strike-slip deformation within the Longitudinal Valley (Chang et al., 2000;Bruce et al., 2006;Shyu et al., 2006;Thomas et al., 2014).The earthquake ground motions were recorded by six near-fault strong motion sensors; station F073 recorded velocities >1 m/s (Yagi et al., 2023) (Fig. 1). Lee et al. (2023) and Tang et al. (2023a) produced kinematic source inversions for the mainshock using teleseismic, strong-motion, and 1 Hz GPS data, both suggesting that the rupture propagated northeastward from the epicenter at ~2.5 km/s on a west-dipping fault.Their models together with aftershock seismicity distributions (Sun et al., 2024) (Fig. S1) suggest however that the east-dipping Longitudinal Valley Fault (LVF) also participated in the earthquake, with localised asperities (<10 km) reaching slip magnitudes of ~1 m.Several bridges close to the earthquake source were damaged (Ko et al., 2023).Two bridges suffered complete structural collapse, which has been attributed to long-period ground motions recorded at nearby strong-motion stations (Carey et al., 2023;Ko et al., 2023) (Fig. 1).
In contrast to previous studies that focused on the broad aspects of the Chihshang earthquake, we focus on the near-fault ground velocities captured by three strong-motion and six 1-Hz GPS sensors (Fig. 1).In particular, we take advantage of a near-fault pair of sensors (250-3800 m from the fault) located on either side of the fault in the northern part of the rupture where surface slip is >1.5 m (we refer to this as the main asperity; black dashed line in Fig. 1).Both of these sensors record dynamic particle displacements during the passage of the rupture front.Using forward dynamic rupture models, we explore the dynamics of the Chihshang earthquake, and show how spatial variations in the shallow oblique slip direction leads to marked asymmetries in near-fault ground motion either side of the fault rupture.

Strong-motion data
In order to derive velocity and displacement time-series data from the strong-motion accelerograms, we first remove the instrument response (instruments contain SMART24A sensors).Next we integrate the acceleration data in time to obtain velocity time-series.Figure 2A shows velocity waveforms for station G020.After analyzing the velocity waveform data closely, we apply a linear correction to the beginning and end of the trace to remove the instrument tilt (Fig. 2A), and then integrate the velocity in time to obtain displacement time-series (Figs.2B and 2C).We use 1 Hz GPS sensor TTCS located <100 m from G020 to compare the magnitude of static displacement (Fig. 2B), and the dynamic particle motion (Fig. 2C).Both aspects of the displacement time-series agree well with the GPS data at this location, giving us confidence that the displacement time-series results are reliable.

Co-seismic displacements from optical image correlation
We use subpixel correlation (Leprince et al., 2007) of Sentinel-2 optical images acquired before (September 12th and 17th, 2022) and after (September 27th and October 10th, 2022) the mainshock to map the horizontal ground displacement during the Chihshang earthquake and to reveal the orientation and extent of the source fault rupture (shading in Fig. 1B).To produce the displacement map, we apply a phase correlation scheme to all image pairs with a correlation window size of 32×32 pixels and a step size of 8 pixels.We obtain a set of displacement maps resolved at 80 m that are then stacked to reduce noise (Fig. 1).The north-south component of the horizontal displacement field reveals a clear and continuous 45 km-long, north-east trending fault rupture within the longitudinal valley, implying slip close to (or intersecting) the ground surface.Displacements reach >1.5 m in the northern part of the rupture near the township of Yuli, where field investigations confirm surface fault rupture and uplift of the western side of the fault (Lee et al., 2023;Yagi et al., 2023) (Fig. S2).
Static displacements from GPS sensors and strongmotion seismometers show subsidence of the Eastern Coastal Range to the east of the fault (dark grey symbols in Fig. 1B), and uplift of the Central Range to the west (white symbols in Fig. 1B) suggesting oblique slip with left-lateral and reverse components on a west-dipping fault, possibly the Central Range Fault.Horizontal displacements from stations on the hangingwall are largely strike-parallel, while stations to the east on the footwall show a dominant fault-normal offset (Fig. 1B).This asymmetry is exemplified by a pair of sensors occupying both sides of the fault rupture above the main asperity-strong-motion sensor F073, located 250 m to the west of the fault, and GPS sensor JPIN, located 3800 to the east of the fault.

Dynamic displacements at F073 and JPIN
Stations F073 and JPIN captured a remarkable nearfield record of dynamic particle motion on both sides of the fault rupture (red and blue curves, Fig. 3A).Strongmotion station F073 (100 Hz) records initial strikenormal displacement towards the fault, followed by a pulse of fault-parallel displacement ~1 m in amplitude, and ending with strike-normal motion in the opposite direction away from the fault.On the opposite side of the fault, particle motion at 1 Hz GPS receiver JPIN is characterized by initial strike-normal displacement towards the fault (~0.3 m), followed by motion in a direction ~45°angle to the strike of the fault.Following this, JPIN particle motions become polarized in the strikeparallel direction, and ~0.4 m of dynamic displacement is abruptly recovered in a 'whiplash' type of motion before the station arrives at its permanent static offset.
These two records show large differences in their dynamic motion, despite being at near-field locations at similar distances along strike of the fault rupture.The displacement history at F073 contains the transient signature of a fault-normal Rayleigh wave and the permanent fault-parallel displacement expected for near-fault strike-slip particle motion (e.g.Aagaard and Heaton, 2004), while JPIN records no Rayleigh wave, and instead shows large fault-normal permanent offset consistent with reverse slip on a nearby fault (e.g.Oglesby and Day, 2001).

Estimating slip-weakening distance from fault-parallel velocity at F073
An important parameter in fault mechanics and physics-based models of dynamic earthquake rupture is the slip weakening distance (D c ), defined as the amount of slip accrued as fault strength drops from static friction to a residual dynamic friction level (Ida, 1972;Palmer and Rice, 1973).Mikumo et al. (2003) showed that the amount of slip (D c ') at the time of peak slip velocity is a close approximation to D c on the fault.Direct on-fault measurement of displacement during seismic slip on natural fault surfaces is challenging; however, velocity seismograms recorded at locations extremely close to the fault rupture (<4 km) closely resemble the on-fault slip velocity on the nearby fault plane (Mikumo et al., 2003;Fukuyama et al., 2003;Fukuyama and Mikumo, 2007).In partic-ular, the amount of displacement at the time of peak particle velocity recorded at near-fault strong-motion stations is used to derive an approximate value of D c , D c " (Fukuyama and Mikumo, 2007).Near-fault seismograms of this type, although still relatively rare, have been used to estimate the on-fault slip weakening process during earthquake ruptures across a range of magnitudes from Mw 6.2-7.9 (e.g., Mikumo et al., 2003;Fukuyama and Mikumo, 2007;Fukuyama and Suzuki, 2016;Kaneko et al., 2017;Chen et al., 2021;Kaneko and Goto, 2022; Cruz-Atienza and Olsen, 2010).
We apply this method to the velocity record of strong- motion station F073, which is located ~250 m from the fault rupture (Fig. 4B).Particle velocity at F073 is characterized by a fault-parallel velocity pulse with peak amplitude of 1.05 m/s and period of ~2 s.Interestingly, the amplitude of peak fault-normal velocity is roughly half (~0.5 m/s) (Fig. 4A).The fault-parallel velocity pulse at F073 is complicated by higher-frequency signals, possibly from site effects, and therefore contains a double peak (Fig. 4B).We use the first of these two peaks, which corresponds to 0.45 m of displacement, to derive a minimum estimate for D c " of 0.9 m (using a simple assumption of an equal contribution from both sides of the fault).The assumption of symmetrical displacement on either side of the fault is not valid for this oblique-slip event, however we use the value of 0.9 m as derived from the fault parallel component, while acknowledging that this is only an approximation and represents a maximum estimate of D c " (because the displacement should be larger on the hanging wall side than on the footwall side due to effect of the fault geometry).We use this value of D c " to inform our selection of slip-weakening distance for the shallow portion of our dynamic rupture model (D c2 ).
3 Dynamic Rupture Models of the Chihshang Earthquake

Model setup
To model the Chihshang earthquake source, we simulate spontaneous dynamic rupture on a 50 km-long, 10 km-wide dislocation embedded within an elastic halfspace.We use an average strike angle of 201°from our optical image correlation results and assign a prestress rake angle of 25°and dip angle of 60°westward from the USGS moment tensor catalogue.Seismic velocity increases in a step-like fashion, approximately following a local 1D velocity model (e.g., Fig. 3 of Huang et al., 2014), and includes a predominant low-velocity layer at shallow depth (Table S1).We use a linear slip-weakening friction law (Ida, 1972;Palmer and Rice, 1973), in which the friction coefficient linearly decreases from its static value μ s to a dynamic value μ d over a characteristic slip distance D c .We assume that the shallower portion of the fault from 0 to 5 km distance along dip is character- As a result, both the static (τ s ) and dynamic (τ d ) strength as well as the magnitude of initial shear stresses τ 0 linearly increase with depth (τ 0 = 0.55σ).We set the value of frictional cohesion to 1 MPa between the free surface and 2 km down dip distance, and to 0 MPa from 2 km down dip distance to the base of the fault model.Dynamic rupture is nucleated at a down-dip distance of 8 km near the southern edge of the fault, and then propagates spontaneously across the fault surface.The numerical code we use is based upon a spectral element method (Ampuero, 2002;Kaneko et al., 2008), which has been verified through a series of benchmark exercises (Harris et al., 2018).

Comparison with observations
First, we vary the dynamic stress drop on the fault by adjusting the value of μ d while keeping μ s =0.67.We match the moment magnitude of the Chihshang earthquake (Mw 7.0) with a value of μ d =0.44.Next we compare synthetic and observed particle motions for stations F073 and JPIN and find that while our initial model with uniform prestress rake of 25°matches the particle motions at F073, the shape of modelled particle displacement does not agree with the observations at JPIN (Fig. 3A).By varying the prestress rake angle between 25°-60°with a different combination of μ s , μ d , D c1 and D c2 , we show that no fault model with a single uniform rake angle ("single patch") can reproduce the static or dynamic motion at both near-fault stations simultaneously.We infer that the particle motions observed at F073 suggest slip at a shallow rake angle (grey curve, Fig. 3A), while the record at JPIN suggests reverse slip at a steep rake angle (black curve Fig. 3A).
Considering that F073 and JPIN are not the same distance from the surface rupture, and that the pattern of surface uplift suggests the fault dips west underneath F073 and away from JPIN, the particle motions at these two stations are likely the result of their sensitivity to fault slip at different depths on the fault surface.Moreover, because JPIN is located km from the fault rupture, it is likely that this down-dip transition in rake angle is shallow.To test this, we consider several multipatch models where pre-stress rake angles change from strike-slip dominated in the shallow part of the fault (0-2 km down dip distance), to reverse slip on the deeper parts of the fault (2-10 km down dip distance) ("multi patch", Fig. 3B).To better match the observations, we shortened the length of the multi-patch fault model from 50 to 46 km.We varied pre-stress rake angles between 0-25°for the shallow patch and between 50-70°for the deeper patch.The best fit to the particle motions was achieved with a pre-stress rake angle of 0°for the shallow patch and 70°for the deeper fault patch.Despite the shallow fault patch having a pre-stress rake angle of 0°, and therefore no component of dip-slip, dynamic stresses from reverse slip at depth cause the co-seismic slip direction to be oblique near the surface (Kearse and Kaneko, 2020), thereby reducing the alongdip contrast in rake angle in the final slip distribution (Figs.3B and S9).To further demonstrate that a shallow transition is a better fit to the observations, we ran the same model but with a deeper (5 km down dip distance) transition in rake angle.This model resulted in a significantly poorer fit to particle motion at JPIN than our preferred model (grey curve in Fig. 3B), due to the dominant influence of the shallow strike-slip prestress patch.The evolution of slip rate and maps of final slip distributions for our preferred multi-patch model and the single-patch model are shown in Figure S9.
Kinematic source models (Lee et al., 2023;Tang et al., 2023a) and aftershock distributions (Fig. S1) of the Chihshang earthquake sequence suggest that slip on the northern part of the LVF near JPIN may have occurred during the mainshock rupture.Slip on the LVF may have contributed to the particle motions recorded at JPIN, which is located ~3 km above the east-dipping LVF plane (Fig. S10).We therefore conducted tests to assess the possible bias introduced by ignoring this contribution in our analysis.We conducted dynamic rupture simulations of slip on an isolated, 70°-dipping 10x10 km fault patch, which results in the distribution of slip similar to the kinematic model of Tang et al. (2023a) (i.e., 1 m of slip at ~8 km down dip distance, tapering to near zero slip at the free surface) (Fig. S11).The resulting particle motions at JPIN are small (<15 cm) (Fig. S10) compared to the observed particle motion (>1 m) demonstrating that slip of the LVF has little influence on ground motions at JPIN, and cannot explain the asymmetry in ground motions across the main surface rupture.
We also compare static offsets derived from GPS and strong-motion stations with the final surface displacement field from our dynamic rupture model (Fig. 3C).The multi-patch model reproduces the observed uniform northwest-directed motion of the footwall side of the fault (grey vectors) and the strike-parallel displacements near the fault on the hangingwall (black vectors).Horizontal displacement vectors on the hangingwall show spatial variation with distance normal to the fault, which matches the observations at distances up to 5 km.Further west into the Central Range hanging wall (5-15 km from the fault), modelled displacements show an increasing component of fault-normal convergence.We cannot ground truth this aspect of the model as there are no available geodetic data in that area.

Parameter uncertainties
To understand the sensitivity of the input parameters on the resulting waveform fits, we used the multi-patch model to perform a set of simulations with a different combination of slip-weakening distance, D c2 , and depth extent of the shallower patch.To quantify the waveform misfit between the model and data, we use the net normalized misfit X, following Goto et al. ( 2019) where d and m are observed and computed waveforms, respectively, |d−m| is the root-mean-square of the resid-uals.Both the observed and computed waveforms are lowpass-filtered at 2 s, and the summation is taken over three components and N stations.Since the fracture energy via slip-weakening distance D c strongly influences the overall rupture speed (e.g., Kaneko and Goto, 2022) we seek to fit the waveform shapes instead of the rupture arrival times.We first calculate the time shift that maximizes the cross-correlation coefficient between observed and synthetic velocity waveforms, and then compute the waveform misfit.Instead of focusing on misfit reduction among many stations (N ≥5), which leads to a strong trade-off between the value of D c2 and its along-dip extent (Fig. S3), we choose to leverage the pair of near-fault stations, F073 and JPIN, to better constrain D c2 at a single location within the main asperity where the sensitivity to shallow D c2 is largest.
Because the two ground motion records at F073 and JPIN are dominated by long-period pulse-like velocities they are likely to be closely related to the slip rate function on the fault at shallow depth, and thus to D c2 .In addition, because our model fault geometry, frictional parameters, and velocity structure are uniform along strike, using data from a narrow along-strike window (<1 km) would reduce the influence on waveform misfit of any unmodelled along-strike heterogeneity in real fault properties between stations.
Although a trade-off between the value of D c2 and its along-dip extent is still evident, we identify a preferred model with the smallest misfit (open circle in Fig. 4C).The overall waveform fits in the preferred model are satisfactory (Fig. 4A), given the relatively simple parameterization of the dynamic rupture model.Remarkably, the preferred model reproduces the phase arrival times well at 5 near-fault stations, and no arbitrary time shifts are applied to the waveforms shown in Fig. 4A.We do not assess the model fits to stations at larger fault distances (>5 km) due to unmodelled 3D velocity effects, although the velocity waveforms for all stations are shown in Figures S5 and S6.The preferred model has D c2 = 0.7 m and its along-dip extent of 5 km, which might signify a significant change in the width of the damaged zone (e.g., Cochran et al., 2009).With fixed D c2 and its depth extent as in the preferred model, we additionally vary non-dimensional prestress τ 0 = τ 0 − τ d / τ s − τ d (Kaneko and Lapusta, 2010) and fault frictional cohesion (Figs.S4a and S4b).While a trade-off exists in these model parameters varied here, we found that the convexity of near-fault particle motions (Fig. 3C) is insensitive to these parameters and hence a robust signature for the spatial variation of rake angles.

Along-dip changes in rake angle in oblique-slip ruptures
Near-fault geodetic measurements of static and dynamic motion during surface-breaking earthquakes provide rare insight into co-seismic slip on the shallow fault surface at high spatial resolution that is not possible with sparse instrumentation.In the case of the 2022 Chihshang earthquake, our modelled particle motions for F073 and JPIN suggest a change in the rake angle of slip on the fault at a depth of ~2 km.This rotation of rake angle probably occurs over some finite width, however we do not explore this in our simple models.It is common for oblique-slip faults in transpressional settings to produce a pattern of co-seismic surface deformation similar to the 2022 Chihshang earthquake.For example, the 2010 Mw 7.0 and 2021 Mw 7.2 earthquakes in Haiti (Raimbault et al., 2023;Hayes et al., 2010), and the Mw 7.0 2022 Abra earthquake in the Luzon Arc, Philippines (Tang et al., 2023b) each produced uplift and fault-parallel displacement of the hangingwall, while accommodating fault-normal displacement of the footwall.The slip distributions for these events consistently show oblique slip rake angles at depth, while showing increasing strike-slip on the fault near the free surface.Although this appears to be a general pattern in oblique convergent settings, we resolve unusually large changes in rake angle on the shallow fault plane in the 2022 Chihshang earthquake.Unlike the events dis-cussed above, detailed distribution of shallow slip during the Chihshang earthquake was able to be resolved due to the near-fault geodetic records of dynamic particle motion, combined with models of dynamic earthquake rupture.
Geological data from oblique-reverse faults suggest that this co-seismic slip behavior can persist over multiple earthquake cycles.Up-dip shallowing of rake angles is observed in geological data from the Clarence Fault, which accommodates oblique convergence in New Zealand's Marlborough Fault System (Nicol and Van Dissen, 2002).The authors attribute this style of long-term slip behavior to a listric fault geometry (steeper fault dip angles near the surface compared to dip angles at depth).In their model, oblique convergence is accommodated at depth by oblique-reverse slip on a moderately dipping fault plane.Because the steeper part of the fault near the surface is more favorably oriented for strike-slip, convergence is accommodated off-fault as distributed shortening and uplift within the hangingwall, while strike-slip occurs on the fault plane.Although our models of the 2022 Chihshang earthquake do not consider non-planar fault geometry, it is possible that the fault surface is curved or kinked and steepens near the surface.Stress orientations computed from seismicity within the Chihshang area between 1991-2007 provide supporting evidence (Wu et al., 2010).Wu et al. (2010) show that between 0-10 km depth, SH is oriented ~30-45°anticlockwise to the strike of the Chihshang earthquake rupture.At depths of 10-30 km, SH orientations are close ~90°to the Chihshang earthquake rupture.This type of depthdependent stress orientation favors reverse slip on a low-dipping fault at seismogenic depth, and oblique slip on a steeper fault surface at shallow depth, similar to our model of the Chihshang earthquake.As this behavior is seen both in recent earthquake ruptures and in long-term geological datasets, this may represent a common way transpression is accomplished on a single oblique-slip fault, rather than by partitioning of strikeslip and shortening on separate structures.

Asymmetry in ground motions across the
Chihshang earthquake rupture Strong differences in near-fault ground motions are known to occur across shallow-dipping thrust fault ruptures that break the free surface (Brune, 1996;Oglesby et al., 1998).Most famously, this effect was observed during the 1999 Mw 7.6 Chi-Chi earthquake which ruptured a low-angle thrust fault (30°dip) and produced larger particle displacements and velocities at the free surface above the hangingwall, relative to the those on the other side of the fault (Oglesby and Day, 2001).Our model of the near-fault rupture dynamics constrained by the data at F073 and JPIN, suggests a different type of ground motion asymmetry may occur across steeply dipping oblique ruptures.Figure 5 depicts the asymmetric ground motions in the vicinity of the propagating rupture in our preferred model.On the hangingwall side of the fault, horizontal ground motions are dominated by large fault-parallel velocities, while the footwall side of the fault experiences largely fault-normal velocities.
Such asymmetry across oblique ruptures has implications for structures built in proximity to active oblique slip faults.For example, during the 2022 Chihshang earthquake, the Gaoliao bridge, located ~4 km north of F073 and JPIN, suffered complete structural collapse during the earthquake due to shear failure of its reinforced concrete piers (Fig 1).Geotechnical reports (Carey et al., 2023;Ko et al., 2023) concluded that strong fault-parallel ground motions recorded nearby at F073 (hangingwall side) were likely the cause for the bridge collapse.Importantly however, the bridge was located entirely on the footwall side of the fault rupture (Fig. 6), and based on our modelling, it probably would not have experienced large amplitude ~1 m/s fault-parallel ground velocity characteristic of hangingwall near-fault ground motions.The bridge may have collapsed due to only ~0.4 m/s of fault-parallel ground velocity as suggested by the simulation (Fig. 6).Based on our models of ground motion due to slip on the Longitudinal Valley Fault (Figs.S10 and S11), we suggest that the LVF slip did not contribute much to the damage of the Gaoliao bridge.The Luntian bridge collapse (Fig. 1) was located 20 km along strike to the southwest, and mainly occupied the hangingwall side of the fault, although the eastern abutment was cut by fault rupture.GNSS station DCHU is located immediately east of the surface rupture near the bridge on-ramp, and therefore recorded only footwall motion.Large differences in ground motion across small distances are known to result from local site effects; however, this work demonstrates how earthquake source dynamics alone can lead to marked asymmetry in strong ground motion at near-fault locations.

Fault-parallel ground velocity pulse
Earthquake ruptures that propagate at velocities above the shear wave speed of crustal rocks (typically >3 km/s), so called 'supershear' earthquakes, are known to produce ground motions that are dominated by a large fault-parallel velocity pulse at near-fault locations, such as the "Pump Station 10" record for the 2002 strikeslip Denali earthquake (Dunham and Archuleta, 2004;Ellsworth et al., 2004).Peak horizontal ground velocity observed in the fault-parallel direction (rather than fault-normal direction typical of rupture directivity) has also been reproduced in supershear laboratory earthquakes (Xia et al., 2004), and has been referred to as a signature of supershear rupture propagation (Mello et al., 2010).Yet, there is mounting evidence that large, pulse-like fault-parallel velocity in excess of fault-normal velocity does occur in the nearfield of earthquake ruptures that do not exceed the shear wave speed (Kaneko and Goto, 2022).Although such observations are rare, ground velocity records of this type have been observed in the 2000 Mw 6.6 Tottori earthquake, Japan (Fukuyama and Mikumo, 2007), 2016 Mw 7.1 Kumamoto earthquake, Japan (Kaneko and Goto, 2022), and the 2016 Mw 7.8 Kaikoura earthquake, New Zealand (Kaneko et al., 2017).Ground velocity recorded at strong-motion station F073 during the 2022 Chihshang earthquake represents another example of this (Fig. 4B).
A near-fault, fault-parallel velocity pulse is interpreted as the result of a fling step effect (Abrahamson, 2001;Kalkan and Kunnath, 2006) where near-fault displacement for a simple, Haskell source model becomes a ramp function with a rise time and the corresponding ground velocity becomes a box car function (Haskell, 1969), resulting in a pulse-like ground velocity.While the fault-parallel velocity pulse is dominated by the static effect at near-field locations, the factors controlling the amplitude and period of the fling step remain unclear.Our models show that at the free surface, large on-fault slip rates (>1 m/s) are sustained for longer time periods (~1 s), compared with slip at greater depth (<0.25 s at 6 km down dip distance) (Fig. 7).We suggest large slip rates are maintained by the dynamic interaction of the propagating rupture and free surface which is enhanced by reflected shear waves from the boundaries of low-velocity layers, as shown by Kaneko and Goto (2022).In addition, shallow rupture speeds in excess of the low velocity material may also play a role in the production of enhanced fault-parallel ground motions (weak supershear signatures can be seen in Figs. 5 and 6).In any case, we suggest that the period and amplitude of near-fault ground motions such as at F073 are controlled by the dynamics of rupture propagation at shallow depth.More investigations are needed to fully untangle the various mechanisms that can en-hance long-period near-fault ground motions.

Constraints on slip-weakening distance
Strong-motion records at near-fault locations during large magnitude earthquakes are still rare, meaning it is important to document each case as it occurs.Our preferred model with the lowest waveform misfit has a D c2 value of 0.7 m, close to D c " = 0.9 m derived from the method of Mikumo et al. (2003) using the record of F073 (Section 2.4).However, there is a clear tradeoff between the width of D c2 and D c2 length, making it difficult to constrain precisely.Moreover, due to the near-fault distance of F073, these results are only relevant for the shallow part of the fault.Our estimate is similar to values of D c " obtained from near-fault highrate GPS records, which show that D c " is ~30-40% of total surface slip.We note that for simplicity, our dynamic rupture model assumes a slip-weakening friction law without accounting for velocity-strengthening friction at shallow depths (Kaneko et al., 2008) or off-fault plasticity (e.g., Andrews, 2005;Ma, 2008;Kaneko and Fialko, 2011).Examining how estimated D c from simplified models, such as the one used in this study, may be influenced by these effects remains to be investigated.

Slip on the Longitudinal Valley Fault during the Chihshang earthquake
Our earthquake rupture model reproduces the static and dynamic geodetic pattern of deformation within Longitudinal Valley without including the secondary, east-dipping Longitudinal Valley Fault (LVF).However, there are reports of ground surface rupture (up to ~20 cm offsets) across the mapped trace of the LVF at the time of the mainshock rupture (Ko et al., 2023).Lee et al. (2023) modelled the Chihshang mainshock earthquake as a west-dipping fault source (similar to our model) but with simultaneous rupture of the east-dipping LVF that contributed 17% of the total seismic moment of this event.Long-period waves recorded at station HGSD, located to the northeast and beyond the extent of our imaged fault source (Fig. 1) were used to argue for significant slip on the LVF (Lee et al., 2023), yet our source model results in ground motions there that match those recorded by HGSD (Fig. S5), suggesting slip on the LVF may not be required to explain the HGSD waveforms.
Although it is likely that slip of <1 m did occur on the LVF, the near-fault strong-motion and geodetic data that we have analysed do not show clear signals from the LVF slip.

Figure 1
Figure 1 Setting of 2022/09/18 Chihshang earthquake.A) Location and extent of main figure within Taiwan.B) Static displacements caused by Chihshang earthquake mainshock from strong-motion accelerometers (squares), 1Hz GPS stations (circles), and north-south component measured from optical image correlation of Sentinel 2 images (shading).Vertical displacements recorded by instruments are indicated by either white (up) or grey (down).Focal mechanisms for the foreshock and mainshock are from the Real-Time Moment Tensor Monitoring System (Lee et al., 2013).Cyan vectors show displacements at co-located GPS TTCS site and strong-motion G020 site (see Figs. 2A-2C).Note pair of sensors F073 (250 m from fault) and JPIN (3800 m from fault) near fault rupture at latitude 23°20'N.Dashed black line near F073 shows location of main asperity where surface slip is >1.5 m.Photographs of collapsed bridges from Ko et al. (2023).

Figure 2
Figure 2 Static and dynamic particle motions from strong-motion seismometers.A) North component of velocity timeseries at strong-motion station G020 (red).Integration to displacement is achieved by first applying a baseline drift correction (black).B) Comparison of displacement time-series with GPS station TTCS (blue) shows that the applied correction is robust.C) Comparison of displacement time-series with TTCS for shorter periods shows that the high velocity (≤0.5 m/s) dynamic particle motion is reliably resolved by this method.
ized by a larger value, D c2 = 0.7 m (similar to our analysis of D c " at F073), while D c1 = 0.3 m on the deeper part of the fault from 5 to 10 km.Model prestress magnitude and rake angle are uniform along strike and prestress magnitude increases linearly with depth.Effective normal stress σ increases with depth via σ = ρgz sin δ(1 − λ) = 10z (MPa), where δ is fault dip angle in degrees, g is gravity in m/s 2 , density ρ = 2700 kg/m 3 , z = distance along dip in km, and λ = 0.62 is the fluid pressure ratio.

Figure 3
Figure3A) Dynamic particle motions for F073 and JPIN compared with modeling results using a single fault patch.In each of the 3 single-patch models displayed, the prestress rake angle is uniform 25°(light gray), 45°(dark gray), and 60°(black).Particle motions for both F073 and JPIN are not well reproduced with any of these models.B) Comparison with a multi-patch model where prestress rake is strike-slip between the free surface and 2 km down dip.In the epicentral region, prestress rake angle is 25°(USGS slip model) at depth, and 70°in the northern part of the fault beneath F073 and JPIN.The downdip variation in rake angle is needed to satisfy particle motion records at both F073 and JPIN.C) Map view of model static displacement field.Observed static offsets (red vectors) shown for comparison.Hanging wall displacements (black vectors) show fault-normal motion away from the fault, and fault-parallel motion near the fault trace.

Figure 4 A
Figure 4 A) Comparison of observed (black) and predicted (red) velocity waveforms calculated with the preferred model for selected stations, both low-pass filtered at 2 second periods.Unfiltered velocity is shown for F073 (magenta).B) D c " estimated from the fault-parallel velocity at F073 from the method of (Mikumo et al., 2003).Fault-parallel displacement is 0.45 m at time of peak velocity.C) Parameter testing results shown as trade-off between magnitude of D c2 in the shallow part of the fault (x axis), and the width of the shallow part of the fault (y axis).Colors represent velocity waveform misfit between the model output and the stations as shown in (A).Preferred model parameters shown with black circle.

Figure 5
Figure 5 Asymmetry in ground motion across the fault in the Chihshang earthquake.Left) vectors show the direction and amplitude (colors) of particle velocity at time steps 15-18s after the origin time.Horizontal line at y=30 km is the fault trace.Red and black sections represent broken and unbroken parts of the fault, respectively, separated by the rupture tip.Data from F073 and JPIN plotted for comparison.Right) particle motions for the same time steps at F073 and JPIN colored by particle velocity.

Figure 6
Figure 6 Snapshot of modelled fault parallel ground velocity at 17.5 seconds.Inset shows the the velocity time history at F073 (grey) and the location of the collapsed Gaoliao bridge (magenta).The model shows large asymmetry in fault-parallel ground velocity across the fault.Lower amplitudes are predicted for the footwall side of the fault where the collapsed Gaoliao bridge is located (magenta) compared to F073 on the other side of the fault.

Figure 7
Figure7Comparison of model slip rate at 6 km down dip distance (black) and model slip rate at the free surface (red) at the location of F073 and JPIN near the main asperity of the Chihshang earthquake.Red curve shows sustained slip at high velocity compared with the more impulsive slip rates at greater depth (black curve).