Multimode ambient noise double-beamforming tomography with a dense linear array: revealing accretionary wedge architecture across Central Taiwan

SUMMARY

Taiwan, one of the most active orogenic belts in the world, undergoes orogenic processes that can be elucidated by the doubly vergent wedge model, explaining the extensive island-wide geological deformation. To provide a clearer depiction of its cross-island orogenic architecture, we apply ambient noise tomography across an east–west linear seismic array in central Taiwan, constructing the first high-resolution 2-D shear velocity model of the upper crust in the region. We observe robust fundamental- and higher-mode Rayleigh waves, with the latter being mainly present in the western Coastal Plain. We develop a multimode double-beamforming method to determine local phase velocities across the array between 2- and 5-s periods. For each location, we jointly invert all available fundamental- and higher-mode phase velocities using a Bayesian-based inversion method to obtain a 1-D model. All 1-D models are then combined to form a final 2-D model from the surface to ∼10 km depth. Our newly developed 2-D model clearly delineates major structural boundaries and fault geometries across central Taiwan, thereby corroborating the previously proposed pro-wedge and retro-wedge models while offering insight into regional seismic hazards.

1 INTRODUCTION

The vigorous orogenic activity in Taiwan is attributed to the oblique convergence between the Eurasian Plate (EP) and Philippine Sea Plate (PSP, Suppe 1984; Teng 1990). This oblique convergence causes intense crustal deformation and mountain building, leading to complex geological structure and multiple tectonic units such as the Coastal Plain, Western Foothill, Hsuehshan Range (HSR), Central Range, Longitudinal Valley and Coastal Range from west to east (Fig. 1). As best depicted by the classic thin-skinned tectonic model, the western part of Taiwan is an accretional wedge formed by a series of folds and thrusts upon a gentle east-dipping decollement (Suppe 1981; Carena et al. 2002; Yue et al. 2005). In contrast, the eastern part of Taiwan accretionary wedge, primarily encompassing the Central Range, undergoes rapid uplift and exhumation. This could involve severe mid-to-lower-crustal thickening through lithospheric collision (Tsai 1986; Yen et al. 1998; Yu & Chen 1994) or underplating (Simoes et al. 2007; Chen et al. 2022). In a context of the doubly vergent orogenic wedge model (Willett et al. 1993; Chuang et al. 2014), the eastern Central Range acts as the retro-wedge in the orogeny system. The eastern boundary of the retro-wedge is named the Central Range Fault (CRF, Biq 1965; Shyu et al. 2006), which plays an important role in mountain building of the Central Range to the west and the formation of the Longitudinal Valley to the east. Based on seismological observations, geodetic measurements and geomorphological mapping, the CRF is suggested to be an active west-dipping reverse and left-lateral fault (Shyu et al. 2006; Wu et al. 2006; Chuang et al. 2014; Huang & Wang 2022), especially after the 2022 September 18 M_w 7.0 Chihshang earthquake in the southern Longitudinal Valley (Lee et al. 2023; Tang et al. 2023; Lin et al. 2024).

The interaction between the EP and PSP has led to complicated structural settings and fault systems in Taiwan, which results in the frequent occurrence of moderate-to-large earthquakes (e.g. Chen et al. 2001b; Yue et al. 2005; Wu et al. 2006; Kuochen et al. 2007; Lee et al. 2014; Huang & Wang 2022). To further investigate the mechanism/evolution of Taiwan's orogenic belt and active fault behaviours, geophysical and geological studies have been conducted including balanced cross-sections (Suppe 1981; Yue et al. 2005), local geological and geomorphological mapping (Johnson et al. 2001; Shih & Teng 2001; Lee et al. 2002), geodetic measurements (Chen et al. 2001a; Hsu et al. 2007; Ching et al. 2011), seismic reflection/refraction imaging (Lai & Chern 2002; Wang et al. 2002; Van Avendonk et al. 2014, 2016; Kuo et al. 2016), geophysical data set analysis (Chang et al. 2000; Kao & Chen 2000), earthquake sequences (Chuang et al. 2014; Huang & Wang 2022), geological trenching/logging (Wu et al. 2007; Hung et al. 2009) and seismic tomography (Kuo‐Chen et al. 2012; Huang et al. 2014a, b, 2015). However, limited depth and lateral resolution of existing regional tomography models prevents robust constraints on the geometry of major faults and structural boundaries. Therefore, constraining detailed subsurface fault geometry remains a major challenge despite its importance to the earthquake hazard assessment.

To better understand the upper crustal structure and orogenic evolution, we apply ambient noise tomography (Shapiro et al. 2005) across an east–west dense linear geophone array in central Taiwan. The array was deployed as part of the TAIGER (Taiwan Integrated Geodynamic Research) project (Fig. 1; Kuo-Chen et al. 2012). Contrary to earthquake-based tomography (e.g. Huang et al. 2014a), ambient noise tomography combined with dense seismic arrays has demonstrated superior resolution on illuminating detailed fault and basin structures at shallow depths (Lin et al. 2013; Wang et al. 2019a; Jia & Clayton 2021; Liu et al. 2021; Rabade et al. 2023). In this study, we first show that robust fundamental- and higher-mode Rayleigh waves can be observed in the noise cross-correlations across the entire central Taiwan linear array and in the western Coastal Plain, respectively. Second, we describe a modified double-beamforming method (Wang et al. 2019a, b; Gkogkas et al. 2021) to resolve multimode Rayleigh wave phase velocities. We then invert Rayleigh wave phase velocity measurements to construct a detailed 2-D shear velocity model across central Taiwan. By comparing the 2-D shear velocity model with the regional seismicity, the new model illuminates subsurface structural geometry and supports the doubly vergent wedge model in Taiwan. The model validates (1) the geometry of the Taiwan Main Detachment beneath pro-wedge in western Taiwan and (2) the existence of the CRF, which intersects with the Longitudinal Valley Fault (LVF) systems and contributes to the uplift of the retro-wedge orogen in eastern Taiwan. The new 2-D shear velocity model also provides new insight into the composition and geometry of major geological structures across the orogen.

2 DATA AND METHODS

2.1 Data pre-processing and ambient noise cross-correlation

We use 73 days of continuous waveform data from 2009 March 20 to May 31, recorded by the linear array of the TAIGER project in central Taiwan (Fig. 1; Okaya et al. 2009; Kuo-Chen et al. 2012). The linear array comprises 56 three-component Mark Products L22 short-period geophones with 2 Hz natural frequency cutting across central Taiwan with station spacing of ∼2 km (Figs 1a and b). We closely follow the data-preprocessing workflow of Liu et al. (2021) on all continuous daily waveforms and construct the multicomponent cross-correlation functions (CCFs). First, we remove instrumental responses, mean and trend from the daily waveforms (Bensen et al. 2007) and decimate the sampling rate to 5 Hz. Then, to preserve the relative amplitude between each component, we filter three-component waveforms between 15- and 50-s periods calculate the absolute mean with a 128-s running time window to obtain temporal normalization functions, and then divide the three-component unfiltered waveforms by the maximum of three-component temporal normalization functions to suppress the undesired energetic earthquake signals (Lin et al. 2014). While 2 Hz geophones are not extremely sensitive to signals above 1-s period, teleseismic signals excited by large earthquakes can still be effectively identified with the 15–50 s bandpass filter (Wang et al 2019b). Next, to whiten the spectrums, we divide the spectrum of each component by the averaged smoothed spectrum of the three-component waveforms. After cross-correlating all the pre-processed daily noise waveforms, we stack and rotate them to obtain nine-component CCFs. In this study, we focus on extracting Rayleigh waves using the vertical–vertical (ZZ) CCFs. In addition, we use the vertical–radial (ZR) component to check the particle motion and exclude non-Rayleigh wave signals. Given that phase velocity measurements exhibit lower uncertainties and possess a deeper sensitivity kernel, revealing deeper structures at the same period compared to group velocity (Lin et al. 2008), this study emphasizes the extraction of phase velocity measurements to construct the first isotropic Vs model across central Taiwan. An example of 2–5 s bandpassed ZZ CCF record section is shown in Fig. 1(c), where both fundamental-mode Rayleigh wave and higher-mode/body wave signals can be clearly observed with waves predominantly propagating eastward (Fig. 1c). In this study, for simplicity, we stack the causal and acausal parts of CCFs to construct symmetrical CCFs for following analyses.

2.2 Iterative double-beamforming tomography

Here, we briefly describe the double-beamforming method developed by Wang et al. (2019b) and only focus on the modification we made to measure both higher- and fundamental-mode Rayleigh wave phase velocities.

We consider all station locations as possible beam centres for both the source and receiver beams. We fix the beam width as 6 km for all the beams to include sufficient stations (> 4), except for one specific beam centred in the East Central Range (black square in Fig. 1b) where a 20 km beam is used due to the station gap. We only keep source–receiver beam pairs with beam centre distance |$\ge $| 3 wavelengths, assuming a reference velocity of 3 km s⁻¹, to satisfy the far-field approximation (Yao et al. 2006). For each source-and-receiver beam pair and each central period, we normalize all available bandpassed CCFs by their maximum amplitudes. We then shift and stack the normalized CCFs, where the time-shift of each CCF is calculated based on the distances between the source and receiver stations to the source and receiver beam centres, respectively, and the assumed source and receiver phase slowness (i.e. inverse of phase velocity). We perform this shift and stack process in the frequency domain, where a time-shift is equivalent to a phase shift, to avoid the finite time-step limitation. We then transform the stacked waveform back to time domain and consider the maximum waveform amplitude as the beam power. We perform a grid search over both source and receiver phase slowness with the slowness ranging between 0.25 and 2 s km⁻¹.

For single-mode process, the double-beam method described by Wang et al. (2019b) would simply take the source and receiver slowness values corresponding to the highest beam power. Here, to account for the multiple energy packages observed (Fig. 1c), we develop an iterative process to separate the fundamental- and higher-mode signals (Fig. 2). In the initial iteration, through grid search, we first determine the source and receiver slowness and the associate group arrival time (⁠|${{{{t}}}_{{\rm{max}}}}$|⁠) that gives the maximum beam power. This corresponds to the primary energy package. We then cut this most energetic phase out from the shifted but unstacked CCFs with a |${{{{t}}}_{{\rm{max}}}}$||$\pm \ 0.15T$| signal window where T represents period. Here, we use the instantaneous period of the stacked CCF based on the phase derivative at |${{{{t}}}_{{\rm{max}}}}$|⁠, which deviates slightly from the centre period of the Gaussian bandpass filter used (Lin et al. 2008). The cut CCFs are then shifted and stacked again in the follow-up iteration to obtain the phase slowness associated with the secondary energy maximum. By comparing the slowness and the |${{{{t}}}_{{\rm{max}}}}$| from the two iterations, we distinguish potential higher- and fundamental-mode energy based on (1) the higher mode should have lower slowness (higher velocity) and (2) the |${{{{t}}}_{{\rm{max}}}}$| of the higher mode should be smaller than the fundamental mode.

We only keep measurements with the stacked CCF signal-to-noise ratio (SNR) above 8, where the SNR is defined as the ratio between the maximum amplitude within the signal window (determined by 0.5 km s⁻¹ minimum and 4.0 km s⁻¹ maximum velocities) and the root-mean-square amplitude of the trailing noise. Furthermore, according to the previous study (Tokimatsu et al. 1992), the fundamental- and higher-mode Rayleigh waves are expected to have a retrograde and prograde particle motion, respectively. Therefore, we implement a selection criterion based on the ellipticity of the particle motion, where only source-and-receiver pairs with a |${{75}^{\rm{o}}} - {{105}^{\rm{o}}}$| phase shift between shift and stacked ZZ and ZR waveforms are kept (Fig. 3). While we observed two energetic phases in most CCFs (Fig. 1c), only stations in the Coastal Plain (green circles in Fig. 1b) show robust higher-mode prograde particle motion. This is consistent with previous studies that observed higher-mode Rayleigh waves in thick sedimentary settings (Tanimoto & Rivera 2005; Savage et al. 2013; Jiang & Denolle 2022). In contrast, stations east of Coastal Plain show linear/tilting instead of prograde/retrograde particle motions, suggesting that the higher-mode energy likely scattered and converted into body waves when entering the mountain terrain (Figs 3c and d). Studying mode conversion is outside of the scope of this study and will be the subject of future investigation.

For each period, beam centre location, and each mode, we estimate the phase velocity and its uncertainty based on the mean and the standard deviation (STD) of the mean of all available measurements (Figs 4b and c). All measurements outside of 2.5 STD are considered outliers and removed from the final statistics. To be conservative and account for unaccounted wave phenomena such as finite-frequency effect, off-great-circle propagation and effect of topography variation, all uncertainties were scaled up empirically by a factor of 1.5 for later Markov Chain Monte Carlo (MCMC) inversion. To highlight the kinks in the observed fundamental-mode phase velocity profiles, for each station location, we calculate the absolute spatial second derivative of phase velocity using the nearest stations on each side (Figs 4e and f). While not exactly one to one, many of the locations with elevated second derivative values coincide with the mapped fault locations and structural boundaries, particularly at the 4-s period.

2.3 Markov Chain Monte Carlo joint inversion

In this study, we follow an MCMC Bayesian inversion method (Shen et al. 2013) to invert available higher- and fundamental-mode dispersion measurements jointly. For each beam centre location, we describe the 1-D Vs velocity model down to a depth of 15 km using six cubic B-splines underlain by a half space model. The Vp and density are set using the empirical relationship described by Brocher (2005). At each location, we extract a 1-D reference model from the 3-D regional model of Huang et al. (2014a) to use as the starting model for the inversion. To stabilize the inversion, in addition to the short-period (2–5 s) dispersion measurements from this study (Fig. 5a), we include 10 to 12-s fundamental-mode Rayleigh wave phase velocities forward calculated using the 1-D reference model. We impose a 10 per cent uncertainty for those ad hoc 10 to 12-s phase velocities, which is larger than the ∼3 per cent uncertainty for the measurements made in this study.

For each location, we randomly perturb the 1-D model (3000 iterations, 16 jumps, up to 75 per cent variation relative to the starting model parameters) and forward-calculate Rayleigh wave dispersion curves (Herrmann 2013). We reject unrealistic models with Vs above 4.9 km s⁻¹ and only accept posterior models with |${{x}^2}$| misfits within |$\pm $| 0.5 of the minimum misfit model. We do not impose any additional weighting between the fundamental- and higher-mode measurements beyond the estimated uncertainties. We then derive the final 1-D Vs model and its uncertainty based on the mean and the STD of the mean of all acceptable models (Figs 7a–c). The predicted phase velocities forward calculated using the inverted model are overall consistent with the measured phase velocities (Fig. 5). We combine the inverted 1-D models from all beam centre locations to construct the final 2-D Vs model.

3 RESULTS

3.1 Phase velocity results

The phase velocity profiles derived from ambient noise cross-correlations and double beamforming (Figs 4 and 5) are generally correlated well with known geological features across central Taiwan. Low phase velocities are observed in the Coastal Plain, particularly at short periods, reflecting the shallow sedimentary deposits. High phase velocities are observed in the Hsueshan Range and the eastern Central Ranges, reflecting the metamorphic core complexes. Slightly moderate fast velocities are observed in the western Central Ranges correlated with where lower grade metamorphic rocks are exposed to the surface. Both the Western Foothills to the west and the Coastal Range to the east exhibit moderate slow velocities corresponding to the sedimentary structure there. These major velocity structures are typically bounded by locations with high second derivative phase velocity values (Fig. 5) suggesting that the velocity changes are sharp. It is worth noting that elevated second derivative values are also commonly observed within the Western Foothills likely reflecting the discrete velocity jumps across the ramp-up faults. We provide a more detailed discussion of the implication of our result in the later section.

It is worth noting that there are significant elevation changes across the Central Range, with a ∼10 per cent inclination on each side (Fig. 4a). In principle, a steep inclination could result in a longer propagation distance for surface waves compared to waves propagating along a flat surface. If the wave is completely bounded to the surface, a 10 per cent inclination would increase the distance by ∼5 per cent, thereby lowering the observed apparent phase velocity by ∼5 per cent. This potential bias is significantly smaller than the > 15 per cent phase velocity variations observed across the Central Range at 2 s (Fig. 4b). Moreover, the observed velocity is lowest where elevation approaches its maximum instead of at locations with the highest inclination. These pieces of evidence suggest that the topography effect is not significant when compared to the intrinsic velocity variation related to subsurface structure. At longer periods, the topography variation (2–3 km) is quite small compared to the wavelength (∼10 km), and hence we expect the potential bias to be even smaller.

To assess the depth sensitivity of our phase velocity measurements, we calculate the 1-D sensitivity kernels numerically by perturbing Vs, Vp and density at different depth (Fig. 6). While it is often assumed that higher-mode Rayleigh waves are more sensitive to deeper structures (e.g. Yoshizawa & Ekström 2010), our sensitivity kernels suggest that it is the fundamental-mode Rayleigh waves that have the deeper sensitivity. Although this is somewhat surprising, it is consistent with the results of Ma et al. (2016), which show that the higher-mode Rayleigh wave eigenfunction could be shallower in a sedimentary setting. This unconventional characteristic likely reflects short-period higher-mode energy being nearly completely trapped in shallow slow sediments due to the sharp velocity increase across the basement. This is consistent with the fact that we do not observe any higher-mode signals outside of the Coastal Plain and beyond 3-s period (Fig. 5).

3.2 Multimode MCMC inversion and model comparison

Fig. 7 shows an example of 1-D inversion with/without including higher-mode dispersion measurements in the Coastal Plain. When using only the fundamental-mode dispersion measurements, a bimodal posterior model distribution between 0 and 3 km depth is observed, demonstrating the trade-off between the velocities at shallow and slightly deeper depths. In this case, the predicted higher-mode dispersion curves (Figs 7b and c) are also bimodal between 2- and 3-s periods, suggesting that the higher-mode dispersion can be used to resolve the tradeoff. By including the higher-mode dispersion measurements, the posterior model distribution is now confined between 0–3 km depth where both the predicted higher- and fundamental-mode dispersion curves (Figs 7b and c) agree with the observations. Benefitting from the complementary sensitivity of higher- and fundamental-mode phase velocity measurements (Fig. 6), the velocity contrast between the shallow slow sedimentary structure and the deeper fast basement can be well resolved across the CP, and the depth resolution could be reliable down to at least 5 km depth. Below 5 km, the broad distribution of posterior models reflects the lack of sensitivity in our short-period measurements (Fig. 6).

The inverted 2-D Vs model (Figs 8a and c) illustrates patterns that are majorly consistent with the observed phase velocity profile (Fig. 5a) and known geological features and tectonic settings (i.e. sedimentary deposits, fault systems and mountain belts). The model also overall agrees with the regional reference model derived from body wave tomography (Fig. 8b; Huang et al. 2014a). Our new model, however, better delineates the geometry of subsurface fine-scale crustal structures. In Figs 8(a) and (c), we superimpose the structural profile of Yue et al. (2005) and the projected background seismicity from 1991 to 2020 onto our model for comparison. Clear correlations can be observed between the inverted velocity structures and known major geological structures including the sharp velocity contrast across the Chinshui Shale beneath the Coastal Plain, the thickening of the foreland fault-and-thrust belt across the Western Foothills, the overall fast metamorphic complex of the Hsuehshan and Central Ranges and shallow sedimentary structure associated with the Longitudinal Valley and Coastal Range.

4 DISCUSSION

While the 2-D Vs model (Figs 8a and c) preserved some of the sharp transitions observed in the phase velocity profile (Fig. 5a), it is less spatially coherent due to minor instabilities in the 1-D MCMC inversions. Consequently, we are unable to identify some of the small-scale velocity transitions in our Vs model using the second derivative analysis, as we did with the phase velocity profile (Fig. 5a). Instead, we present the velocity contours (Fig. 8a) of the spatially smoothed Vs model, with a grid size of approximately 0.023° to represent the average station distance, highlighting the velocity variation. To only focus our discussion on the most reliable part of the model, we masked out depths where the Vs sensitivity of the 5-s phase velocity is less than 0.001 km⁻¹ in our sensitivity kernel analysis (Fig. 6). Generally, the sensitivity kernel tends to move shallower in low-velocity zones due to shorter wavelengths.

4.1 The pro-wedge architecture and pre-orogenic normal faults

Our model further shows the structural details within the doubly vergent accretionary wedge of Taiwan. The western half of the Taiwan island, including Western Foothills, HSR and western part of Central Range, is considered as the pro-wedge, with fold-and-thrust ramping faults propagating above the detachment (Willett & Brandon 2002; Yue et al. 2005; Fuller et al. 2006). While the slow anomaly from our velocity model (Figs 8a and c) shows an excellent geometrical correlation to the geometry of the foreland basin extracted from Lin & Watts (2002) and fold-and-thrust ramping fault systems within Western Foothill from Yue et al. (2005), it also provides more details for understanding the geometry of these structures. We consider Vs 3.0 km s⁻¹ (e.g. ≈ Vp 5.2 km s⁻¹) as a reference velocity separating crystalline from clastic sedimentary rocks (Brown et al. 2022). Background seismicity is mainly observed beneath the Vs 3.0 km s⁻¹ contour line, while the pro-wedge's fold-and-thrust system is primarily developed above this Vs 3.0 km s⁻¹ velocity contour. This correlation reflects the geometry of Taiwan Main Detachment and the thin-skinned thrust model in western Taiwan.

The shallow velocity in the top 3 km in the pro-wedge generally increases eastward, where apparent velocity increases can be observed across the Changhua Fault (CF), Shuangtung Thrust (ST) and the Shuilikeng Fault (SLKF). Layered unconsolidated sediments with velocities slower than 1.0 km s⁻¹ (Fig. 8a) are only observed near the surface west of the CF. This suggests that the CF and its hangingwall anticline served as the deformation front in the pro-wedge, where the foreland basin sediment has been uplifted and eroded (Yue et al. 2005). A similar shallow velocity increase can be found across the ST, where velocity in the top 2 km shows a sharp eastward increase. Note that the velocity change across the ST is much sharper than the change across the CF, suggesting the geometry of the ST is likely steeper than the CF.

In contrast to the shallow velocity increases, step-like eastward velocity drops are observed below ∼4 km depth beneath the CF, ST and Chelungpu Fault (CLPF). These eastward velocity drops have yet to be noticed in the structural profile at this latitude but roughly coincide with where the CF, CLPF and ST are expected to ramp up from the basal detachment. Care must be taken when interpreting these deeper features, as some might reflect trade-off between shallow and deeper structures in the 1-D inversions. At 5-s period, an eastward phase velocity drop can be observed beneath CF in the phase velocity profile (Fig. 5) indicating that at least Vs drop there is a robust feature. We suggest that the divergence of shallow and deep velocity variation hints pre-existing east-dipping normal faults, which were reactivated as reverse faults during the orogeny (Brown et al. 2012). Based on the industrial seismic reflection data, Yang et al. (2016) suggest similar pre-orogenic normal faults exist underneath the Coastal Plain. Simoes et al. (2007) also suggest a high-angle hinge fault exists below the northern part of CF and its anticline, showing similar geometry as we see in our velocity model across the southern part of the anticline. Our finding, along with other lines of evidence from previous analyses, suggests that the pre-existing normal faults could play a significant role in controlling the structural location and geometry of Taiwan's pro-wedge.

Further east, a sharp velocity increase is observed from the surface down to ∼10 km across the surface-mapped SLKF, which is a tectonostratigraphic feature separating the sedimentary Western Foothills and metamorphized HSR. The trend of velocity increase across the SLKF, however, is different from the east-dipping trend across the other faults in the frontal part of the pro-wedge, suggesting the SLKF might be west dipping in contrast to those aforementioned structures. Even though previous studies agreed that the SLKF is the primary feature that separates the synorogenic foreland basin from the metamorphosed zone, the kinematics and detailed structure of the SLKF are still debated and have been explained in various interpretations. For instance, Yue et al. (2005) described the SLKF as a pre-existing west-dipping extensional fault; Brown et al. (2012) exerted it as a steeply east-dipping transpressional fault that penetrates well into the middle crust; Camanni et al. (2014) suggest it as an east-dipping rifted-related fault. From our model, we found that the geometry of the SLKF is consistent with the model proposed by Yue et al. (2005) that the SLKF is a high-angle west-dipping fault and could provide structural evidence of pre-rifting tectonic evolution in Taiwan.

4.2 The retro-wedge architecture

On the eastern side of Taiwan, the eastern part of the Central Range and the Longitudinal Valley are categorized as the retro-wedge of the Taiwan orogen, where the reverse CRF acts as the back thrust structure of the wedge (e.g. Chuang et al. 2014). On the easternmost part of the 2-D model, we observe a distinct velocity contrast between the Central Range and the Coastal Range, where the high-velocity Central Range is bounded by a steeply west-dipping boundary beneath the range. This velocity contrast matches well with the location of the Longitudinal Valley Suture, which consists of the active high-angle west-dipping CRF and the east-dipping LVF (Shyu et al. 2005, 2006). The background seismicity shows a west-dipping cluster and an east-dipping cluster, corresponding to the velocity contrast and delineating the CRF to the west and the LVF to the east. We observed a junction formed between the LVF and CRF approximately at 3 km depth along our velocity profile. This observation is consistent with the interpretation that the CRF is buried beneath the LVF (Shyu et al. 2005, 2006; Chuang et al. 2014; Lee et al. 2014) and subsequently truncates the LVF in the northern Longitudinal Valley (Huang & Wang 2022). The truncation of CRF–LVF plays a crucial part in controlling the coseismic fault propagation for earthquakes generated by the LVF. Also, the velocity contours beneath the surface location of the LVF (Shyu et al. 2020) reveal a steep east-dipping velocity contrast, correlating with the imbricate thrusting geometry mentioned in Hsieh et al. (2020). Combining the previous geological and seismological observations (Shyu et al. 2006; Wu et al. 2006; Chuang et al. 2014; Huang & Wang 2022) with our inferred geometry of the CRF provides additional structural details between the retro-wedge and Longitudinal Valley suture.

4.3 Fault-and-fold structures within the Central Range

We observe a high-angle east-dipping velocity contrast across the Tili Thrust (TT) within the metamorphic HSR area. The TT is a major structural boundary at the western margin of the metamorphic Central Range (Fig. 8c). Previous studies suggested that the TT is a steeply dipping fault, where the velocity boundary can be interpreted as a hangingwall ramp against a footwall ramp (Yue et al. 2005; Brown et al. 2012). Our velocity model also suggests that the TT is a crustal-scale geological boundary that extends to at least 10 km below the surface, with its geometry matching the interpreted structure of the balanced structural profile (Yue et al. 2005). Interestingly, the crustal-scale TT is not correlated to any background seismicity from Wu et al. (2008). Instead, a seismicity cluster with a similar geometry to that of the TT and its deeper ramp appears underneath the SLKF, very close to the eastward velocity drop described in the previous section. Although our profile does not extend deep enough to observe the velocity change associated with this group of seismicity, Chuang et al. (2013) suggest this east-dipping seismicity may be associated with the newly developed duplex system beneath the orogeny belt. While our data cannot prove this hypothesis, like other velocity changes that occurred along the reverse faults within the fold-and-thrust belt in WF, we observe that the location where the seismicity joins the basal detachment aligns with the velocity change in our model. This suggests that the east-dipping seismicity may be associated with a major pre-existing structure within the orogenic belt.

It is worth noting that we also observe the first-order structural geometry within the slate-dominated Central Range. The comparison with the structural profile from Yue et al. (2005) indicates that the two relatively low-velocity zones in the western flank of the Central Range are likely associated with synclinoria, where their cores are composed by the Chiayang formation. Also, at ∼25 km east of the TT, the bowl-shaped velocity drop could imply an unconformity between the Lushan and Pilushan Formations where Mesozoic and Paleozoic schist and marble formations are exposed on the surface along the eastern flank of the Central Range.

5 CONCLUSIONS

In this study, by applying ambient noise double-beamforming tomography across a dense linear geophone array, we construct a detailed 2-D regional crustal shear wave velocity model across central Taiwan. We extract robust Rayleigh wave signals from multicomponent noise cross-correlations, propose an iterative double-beamforming scheme to measure higher- and fundamental-mode phase velocities across the array, and then apply MCMC inversions to construct the final Vs model. The shallow sensitivity of our measurements unveils detailed upper crustal structures with unprecedented resolution. The new model overall is consistent with prior geological and geophysical studies but provides further constraints on the geometry of sedimentary foreland basin and fault systems in central Taiwan. Moreover, the new model provides a clear view of the doubly vergent wedge architecture of Taiwan. For the pro-wedge to the west, the Vs structure reveals detailed geometry of the fold-and-thrust belt and the pre-existing normal faults that may have been reactivated as reverse faults during the orogeny. For the retro-wedge to the east, the CRF is identified as a west-dipping slow anomaly, which is truncated by the east-dipping LVF at ∼3 km depth. These new findings, combined with previous geological and geophysical observations, illuminate the orogenic processes in Taiwan and provide insight into regional earthquake hazards.

Acknowledgement

This work was supported by the National Science Foundation of the US (grant EAR-1753362) and the National Science and Technology Council of Taiwan (grant NSTC111-2116-M-001–032). We thank the TAIGER working group for collecting the data, and John Suppe and Yu-Huan Hsieh for valuable discussion about Taiwan's tectonic and geological structures. We also thank editor Gabi Laske and two anonymous reviewers for their helpful comments.

DATA AVAILABILITY

TAIGER data are accessible on the International Federation of Digital Seismograph Networks (https://doi.org/10.7914/SN/2F_2009). The code to perform double beamforming is available on GitHub (https://github.com/andrewntu/Double-Beamforming.git).

References